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ABSTRACT 

', We present a semi-analytic model to investigate the merger history, destruction rate, and survival 

(•~^ ■ probability of substructure in hierarchically formed dark matter halos, and use it to study the substruc- 

I ture content of halos as a function of input primordial power spectrum. For a standard cold dark matter 

CsJ . "concordance" cosmology (ACDM; n = 1, as = 0.95) we successfully reproduce the subhalo velocity func- 

,— H ' tion and radial distribution profile seen in N-body simulations, and determine that the rate of merging 

^ I and disruption peaks ^10—12 Gyr in the past for Milky Way-like halos, while surviving substructures 

■ are typically accreted within the last ~ — 8 Gyr. We explore power spectra with normalizations and 
' spectral "tilts" spanning the ranges erg ~ 1 — 0.65 and n ~ 1 — 0.8, and include a "running-index" 
, model with dn/dlnfc = —0.03 similar to the best-fit model discussed in the first-year WMAP report. 

We investigate spectra with truncated small-scale power, including a broken-scale inflation model and 

I three warm dark matter cases with mw = 0.75 — 3.0 keV. 

, We find that the mass fraction in substructure is relatively insensitive to the tilt and overall normaliza- 

Q,^ ' tion of the primordial power spectrum. All of the CDM-type models yield projected substructure mass 

I fractions that are consistent with, but on the low side, of published estimates from strong lens systems: 

■ /g = 0.4 — 1.5% (64 percentile), for subhalos < 10^ Mq within projected cylinders of radius r < 10 
' kpc. Truncated models produce significantly smaller fractions, /g = 0.02 — 0.2% for mw — 1 keV, and 

CO . are disfavored by lensing estimates of substructure. This suggests that lensing and similar probes can 

' provide a robust test of the CDM paradigm and a powerful constraint on broken-scale inflation/ warm 
particle masses, including masses larger than the ~ 1 keV upper limits of previous studies. We compare 

Q^. our predicted subhalo velocity functions to the dwarf satellite population of the Milky Way. Assuming 

I ' dwarfs have isotropic velocity dispersions, we find that the standard n — I model over-predicts the 

2 . number of Milky Way satellites at Vmax~ 35 km s~^, as expected. Models with less small-scale power 

^ ' do better because subhalos are less concentrated and the mapping between observed velocity dispersion 

^ I and halo Vmax is significantly altered. The running-index model, or a fixed tilt with cg ~ 0.75, can 

' account for the local dwarfs without the need for differential feedback (for Vmax^ 20 km s^^); however, 
these comparisons depend sensitively on the assumption of isotropic velocities in satellite galaxies. 
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1. INTRODUCTION & Yang 2003), and thc Local Group dwarf galaxy count is 

T,i,ij 1-1 jirxj. r significantly below what might naively be expected from 

in the standard cosmological model oi structure torma- , , , , , . . r ^ ^t^, ^ / i /t/i ■ . i 

i- /A/^T^A4^^ ii, TT • -J • J- J u ij u- ths suDstructurc content ot ACUM haios (Klypm et al. 
tion (ACDM), the Universe is dominated by cold, colli- , t,^ ^ , ^nnn \ t ri ^ o 

■ / J 1 /r^T->i\^\ A a u 1 • 1999a, K99 hereafter; Moore et al. 1999a). In Zentner & 

sionless dark matter (CDM), made Hat by a cosmologi- „ „ , l, r^n^^N i i ,i ; ;i 

, 4. 1. / \\ A A A -t-i, ■ ^ A -J- Bullock (2002, hereafter ZB02), we showed that the cen- 

cal constant (A), and endowed with initial density per- , , , ' . , , , , , , , , 

• , a . .■ A ■ ■ a tral densities ot ACDM dark matter halos can be brought 

turbations via quantum fluctuations during inflation. Ihe . ^ .,, . r i i 

AriTMV/T J 1 -ii, o 1 r> no u n -7 a mto reasonable agreement With the rotation Curves ot dark 
ACDM model with iIm = 1 — iU = 0.3, h « 0.7, and a , • , , , • , , • , ,• , n 

,• ., , r ■ 1-1 i-ui- inti\ matter-dominated galaxies by reducing galactic-scale fluc- 
scale-mvariant spectrum ot primordial perturbations (i-^(fc) (X . . , . °. , •' , / rv r,r ^ 

,„ , r, n^ • 1 ui f 1 i A tuations m the initial power spectrum (ug ~ 0.75 and 

A: , n = 1, cTg ~ 0.9) IS remarkably successful at reproduc- . , ^, , n ^ d txr • 

1 ; 1 r , 1 u / c 1 i n '--^ 0.9 IS a good match; see Alam, Bullock, & Wem- 

mg a plethora ot large-scale observations (e.g., Spergel et , r^^^r, i n. Amur iv«^ ^ , . , ^^^o , 

1 onoQ 13 • 1 f 1 onno^ t f f in berg 2002, hereafter ABW; McGaugh et al. 2003; van den 

al. 2003; Percival et al. 2002 . In contrast, several small- i , i c^r^r^o\ mi , ■ . ■ r 

,1 ,. 1 A-ca ^^ ^ 1 • Bosch et al. 2003 . ihe present paper is an extension ot 

scale observations have proven more difficult to explain. . , ' , , ^ ^ . . ... , 

^ , J ., . J i i 1. 1 this work. We explore how changes m the initial power 

Galaxy densities and concentrations appear to be much , rr-.ii i. . ..r* ^t^h t f ^ 

, ,1 1, i • J- 4- J ^ 4-1, i J J / iN spectrum attect the substructure content ot ACDM halos, 
lower than what is predicted lor the standard (n = l) ^ „ ,. ., ii; 

ACDM model {e.g., Debattista & Sellwood 2000; Cote, ^''f'^'f ^S^?'^^^ attempts to measure the substruc- 

Carignan, & Freeman 2000; Borriello & Salucci 2000; Bin- t^rc mass fraction via gravitational lensing, and relate our 

p t:^ onm j- r,nni J -n> u c rcsults to tflc qucstioii ot thc abuudaucc OI dwait satellites 

ney & Evans 2001; Keeton 2001; van den Bosch & Swa- Local Grou 

ters 2001; Marchesini et al. 2002; Swaters et al. 2003; oca roup. „„, , , , 

\K 1, TD 1 i> J TD1 1 onno A td i, ai^ it is straightforward to see why CDM halos are expected 

McGaugh, Barker, & de Blok 2003; van den Bosch, Mo, , , , ; ; , , n i- ■• . i i i ; 

to play host to a large number oi distinct, bound substruc- 

Hubble Fellow tures. Or "subhalos." in the modern picture of hierarchical 
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structure formation (White & Rees 1978; Blumenthal et al. 
1984; KaufFmann, White, & Guiderdoni 1993) low-mass 
systems collapse early and merge to form larger systems 
over time. Small halos collapse at high redshift, when the 
universe is very dense, so their central densities are cor- 
respondingly high. When these halos merge into larger 
hosts, their high densities allow them to resist the strong 
tidal forces that act to destroy them. While gravitational 
interactions do serve to unbind most of mass associated 
with merged progenitors, a significant fraction of these 
small halos survive as distinct substructure. 

Our understanding of this process has increased dra- 
matically in the last five years thanks to remarkable ad- 
vances in N-body techniques that allow the high force and 
mass resolution necessary to study substructure in detail 
(Ghigna et al. 1998, 2000; Kravtsov 1999; K99; Klypin et 
al. 1999b; Kolatt et al. 1999; Moore et al. 1999a,b; Font 
et al. 2001; Stoehr et al. 2002). For n = 1, ACDM and 
CDM simulations, the total mass fraction bound up in sub- 
structure is measured at / ~ 5 — 15% (Ghigna et al. 1998, 
Klypin et al. 1999b), with a significant portion contributed 
by the most massive subsystems, dN/dM oc M~^, s « 1.7. 
The substructure content of halos seems to be roughly self- 
similar when subhalo mass is scaled by the host halo mass 
(Moore et al. 1999a) and the subhalo count is observed 
to decline at the host halo center, where tidal forces are 
strongest (Ghigna et al. 1998; Colin et al. 2000b; Chen, 
Kravtsov, & Kceton 2003). 

Unfortunately, studies of substructure using N-body sim- 
ulations suffer from issues of numerical resolution. Simula- 
tions with the capability to resolve substructure are com- 
putationally expensive. They cannot be used to study the 
implications of many unknown input parameters and can- 
not attain both the resolution and the statistics needed to 
confront observational data on substructure that appear 
to be on the horizon. Even state of the art simulations 
face difficulties in the centers of halos where "overmerg- 
ing" may be a problem {e.g., Chen ct al. 2003; Klypin et 
al. 1999b) and measurements of the substructure fraction 
via lensing are highly sensitive to these uncertain, central 
regions. Our goal is to present and apply a semi-analytic 
model that suffers from no inherent resolution effects, and 
is based on the processes that were observed to govern sub- 
structure populations in past N-body simulations. This 
kind of model can generate statistically significant pre- 
dictions for a variety of inputs quickly and can be used 
to guide expectations for the next generation of N-body 
simulations. Conversely, this model represents in many 
ways an extrapolation of N-body results into unexplored 
domains and it is imperative that our results be tested by 
future numerical studies. In the present paper, wc aim to 
explore the effect of the power spectrum on the population 
of surviving subhalos, but in principle these methods are 
suitable for testing substructure ramifications for a variety 
of cosmological inputs. 

One of the main motivations for this work comes from 
simulation results that indicate that galaxy-sized CDM 
halos play host to hundreds of subhalos with maximum 
circular velocities in the range 10 km s""^ ^ Knax ^ 30 
km s~^. The Milky Way, as a comparative example, hosts 
only 11 dwarf satellites of similar size. This "dwarf satellite 
problem" specifically refers to the gross mismatch between 



the predicted number of ACDM subhalos and the count of 
satellite galaxies in the Local Group (K99; Moore et al. 
1999a; Font et al. 2001; see also Kauffmann et al. 1993, 
who indicated that there may be a problem using analytic 
arguments). The dwarf satellite problem and other small- 
scale issues led many authors to consider modifications to 
the standard framework. If the dark matter were "warm" 
(Pagels & Primack 1982; Colombi, Dodelson, & Widrow 
1996; Hogan & Dalcanton 2000; Colin et al. 2000a; Bode, 
Ostriker,&Turok 2001; Lin etal. 2001; Knebe et al. 2002) 
or if the primordial power spectrum were sharply trun- 
cated on small scales (Starobinsky 1992; Kamionkowski 
& Liddle 2000) then subgalactic-scale problems may be 
allayed without vitiating the overall success of ACDM on 
large scales. Another possibility is that CDM substructure 
is abundant in all galaxy halos, but that most low-mass 
systems are simply devoid of stars. An intermediate so- 
lution may involve a simple modification of the assumed 
primordial spectrum of density perturbations that gradu- 
ally lowers power on galactic scales relative to the horizon, 
e.g., via tilting the power spectrum. 

Probing models with low galactic-scale power is moti- 
vated not only by the small-scale crises facing standard 
ACDM but also by more direct probes of the power spec- 
trum. While many analyses continue to measure "high" 
values for cts ~ 1 (Van Waerbeke et al. 2002; Komatsu 
& Scljak 2002; Bahcall & Bode 2003; where ag is the lin- 
ear, rms fluctuation amplitude on a length scale of 8 h~^ 
Mpc ), numerous recent studies relying on similar tech- 
niques, advocate rather "low" values of erg ^ 0.7 — 0.8 
(Jarvis et al. 2003; Bahcall et al. 2003; Schuecker et al 
2003; Pierpaoh et al. 2003; Viana et al. 2002; Brown et al. 
2002; Allen ct al. 2002; Hamana ct al. 2002; Mclchiorri 
& Silk 2002; Borgani et al. 2001). Similarly, the Ly-a 
forest measurements of the power spectrum are consistent 
with reduced galactic-scale power (Croft et al. 1998; Mc- 
Donald et al. 2000; Croft et al. 2002). Set against the 
normalization of fluctuations on large scales implied by 
the Cosmic Background Explorer (COBE) measurements 
of cosmic microwave background (CMB) anisotropy (Ben- 
nett et al. 1994), these data suggest that the initial power 
spectrum may be tilted to favor large scales with n < 1. 

The recent analysis of the Wilkinson Microwave Anisotropy 
Probe (WMAP) measurements of CMB anisotropy pre- 
sented by Spergel et al. (2003; see also Verde et al. 2003; 
Peiris ct al. 2003) returns a best-fit spectral index to a pure 
power law primordial spectrum of n = 0.99 ± 0.04 when 
only the WMAP data are considered. However, when data 
from smaller scale CMB experiments, the 2dF Galaxy Red- 
shift Survey, and the Ly-a forest are included, the analysis 
favors a mild tilt, n = 0.96 ± 0.02. Interestingly, all of the 
data sets together yield a better fit if the index is allowed 
to run: the WMAP team find dn/dlnk = -0.03llg:gi?. 
This result is consistent with no running at ^ 2a, and the 
statistical significance is further weakened when additional 
uncertainties in the mean flux decrement in the Ly-a for- 
est are considered (Seljak, McDonald, and Makarov 2003; 
Croft et al. 2002), yet such a model certainly seems worth 
investigating, especially in light of the small-scale difficul- 
ties it may help to alleviate. 

We explicitly show how models with reduced small-scale 
power are expected to help the halo density problem in 
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Fig. 1. — Galaxy central densities. The symbols show the mean 
density (relative to the critical density) within the radius where each 
rotation curve falls to half of its maximum value, inferred from the 
measured rotation curves of several dark matter-dominated dwarf 
and low surface brightness galaxies (see ZB02 for details). The data 
are taken from de Blok, McGaugh, & Rubin (2001; triangles and 
hexagons), de Blok and Bosma (2002; squares), and Swaters (1999; 
pentagons). The lines show the theoretical expectation for several 
of the power spectra we describe in §3 and Table 1. The points with 
error bars in the upper right corner reflect the expected theoretical 
scatter in the density as inferred by Bullock et al. (2001, larger 
range) and Jing (2000, smaller range). 



Figure 1, which is an updated version of Figure 5 in ZB02. 
Here, we compare the densities of standard n = 1 halos 
to galaxy rotation curve data (see ZB02 for details) along 
with expectations for the running-index (RI) model fa- 
vored by WMAP and several other models we explore in 
the following sections. Galaxy and halo densities (vertical 
axis) are evaluated at the radius where the rotation curve 
falls to half its maximum value and expressed in units of 
the critical density {Ay/2, as defined in ABW). Clearly, 
the data favor low small-scale power relative to the stan- 
dard n — 1 case. 

The possibility of discriminating between standard ACDM 
and several alternatives has inspired efforts to measure 
and quantify the substructure content of galactic halos. 
One method relies on studying tidal tails associated with 
known Galactic satellites (Johnston, Spergel, & Haydn 
2002; Ibata et al. 2002a, 2002b; Mayer et al. 2002). 
Subhalos passing through cold tidal streams scatter stars 
away from their original orbits, and the signatures of these 
events may be detectable in the velocity data of future as- 
trometric missions and several deep halo surveys that will 
soon be completed. Of particular interest for obtaining 
measurements in distant galaxies are studies that aim to 
detect substructure via flux ratio anomalies in strong grav- 
itational lenses (Moore et al. 1999a; Metcalf & Madau 
2001; Metcalf & Zhao 2002; Bradac et al. 2002). Us- 



ing a sample of seven four-image radio lenses, Dalai & 
Kochanek (2001, DKOl hereafter) estimated a mass frac- 
tion of / = 0.006 - 0.07 (90% confidence level) bound up 
in substructure less massive than ^ lO'* — 10^" Mq, in 
line with the rough expectations of CDM.^ While mea- 
surements of this kind are susceptible to potential degen- 
eracies with the adopted smooth lens model and other un- 
certainties, they are encouraging and serve as prime mo- 
tivators for this work (see Kochanek & Dalai 2003). In 
addition, new observational techniques that focus on astro- 
metric features (Metcalf 2002), and particularly spectro- 
scopic studies of strong lens systems (Moustakas & Metcalf 
2003), promise to hone in on the masses of the subclumps 
responsible for the lensing signals. 

If the Milky Way really is surrounded by a large num- 
ber of dark subhalos, the dwarf satellite problem serves as 
a conspicuous reminder that feedback plays an important 
role in hierarchical galaxy formation. Of course, the need 
for feedback in small systems has been generally recog- 
nized for some time {e.g., White & Rees 1978). Supernova 
blow-out likely plays a role in regulating star formation 
if CDM is the correct theory (Dekel & Silk 1986; Kauff- 
mann et al. 1993; Cole et al. 1994; Somerville & Primack 
1999); however, supernova winds do not naturally sug- 
gest a sharp discrepancy at ~ 30 km s^ , nor do they ex- 
plain why some halos of this size should have stars while 
most have none at all. It seems more likely that super- 
novae play an important role in setting scaling relations 
in slightly larger galaxies (Vmax ~ 100 km s~^; Dekel & 
Woo 2002; but see Mac Low & Ferrara 1999). Perhaps 
a more natural feedback source on satellite galaxy scales 
is the ionizing background, which should suppress galaxy 
formation in halos with Vmax^ 30 km s^ (Rees 1986; Ef- 
stathiou 1992; Kauffmann et al. 1993; Shapiro, Giroux, & 
Babul 1994; Thoul & Weinberg 1996; Quinn, Katz, & Ef- 
stathiou 1996; Gnedin 2000). Bullock, Kravtsov, & Wein- 
berg (2000, BKW hereafter) suggested that dwarf galaxies 
may be associated with small halos that collapsed before 
the epoch of reionization, and though the method used 
by BKW to estimate dwarf luminosities was crude, more 
sophisticated models have since led to similar conclusions 
(Chiu, Gnedin, & Ostriker 2001; Somerville 2002; Ben- 
son et al. 2002). For the smallest systems, Vmax^ 10 — 
20 km s^^, the ionizing background likely stops star for- 
mation altogether by photo-evaporating gas in halos, even 
after they have collapsed (Barkana & Loeb 1999; Shaviv 
& Dekel 2003). 

Precisely what can be learned about galaxy formation 
and/or cosmology by counting dwarf satellites depends 
upon our expectations for the density profiles of their host 
halos. To count satellites of a given maximum circular 
velocity, we must infer a halo Vmax using the observed 
central velocity dispersion cr*, and the mapping between 
these two quantities depends sensitively on the structure 
of each satellite's dark matter halo. This point was first 
emphasized by S. D. M. White at the Summer 2000 Insti- 
tute for Theoretical Physics Conference on Galaxy Forma- 
tion and Evolution.^ Standard ACDM halos with Knax ~ 

^ DKOl quoted an approximate upper mass limit of 10^ — 10^ Mq . 
They have since concluded that an upper limit of ~ 10* — 10^" Mq 
may be more appropriate (N. Dalai, private communication). 

See http://onlinc.kitp.ucsb.edu/online/galaxy_cOO/whitc/. 
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30 km s^^ are expected to be very concentrated (Colin et 
al. 2000a; Bullock et al. 2001), with their rotation curves 
peaking at rmax ~ 1 kpc, so the multiplicative factor that 
converts a ~ 1 kpc velocity dispersion measurement to 
halo V^ax is fairly modest: Vmax ^ VScTi, (K99). How- 
ever, as we discuss in §4, the appropriate conversion is 
cosmology-dependent because models with later structure 
formation tend to produce halos with more slowly rising 
rotation curves. If a dwarf galaxy sits in a halo with a 
slowly rising rotation curve that peaks at rmax 3> a few 
kpc, the conversion factor, and thus Knax, can be signif- 
icantly larger. Shifts of this kind in the "observed" ve- 
locity function change the implied velocity (or mass) scale 
of discrepancy, and influence our ideas about the type of 
feedback that gives rise to the mismatch. 

Hayashi et al. (2003, H03 hereafter) and Stoehr et al. 
(2002, S02 hereafter) suggested that substructure halos 
experience significant mass redistribution in their centers 
as a result of tidal interactions and that they are therefore 
less concentrated than comparable halos in the field. They 
argue that when this is taken into account, the dwarf satel- 
lite mismatch sets in at Knax ~ 20 km s~^, and that the 
transition is sudden — below this scale all halos are devoid 
of observable galaxies. While these conclusions have yet 
to be confirmed and are dependent upon subhalo merger 
histories and the isotropy of dwarf velocity dispersions, 
they highlight the need to refine our predictions about 
halo substructure. They also motivate us to explore how 
minor changes in cosmological parameters can influence 
our interpretation of the dwarf satellite problem. 

In the remainder of this paper we present our study of 
CDM substructure. In §2, we describe our semi- analytic 
model, provide some illustrative examples, and compare 
our results for standard ACDM to previous N-body re- 
sults. In §3, we briefly describe the input power spectra 
that serve as the basis for this study. In §4, we present 
our results on subhalo mass functions and velocity func- 
tions. We make predictions aimed at niciasuring substruc- 
ture mass fractions via gravitational lensing and address 
the dwarf satellite problem in light of some of our findings 
in this section. In §5 we discuss some shortcomings of our 
model and how they might be improved in future work. 
In §6 we summarize our work and draw conclusions from 
our results. In this study we vary only the power spectrum 
and work within the context of the so-called "concordance" 
cosmological model with flu = 0.3, Ha = 0.7, h = 0.72, 
and ^Bh"^ = 0.02 {e.g., Turner 2002; Spergel et al. 2003). 

2. MODELING HALO SUBSTRUCTURE 

In order to determine the substructure properties of a 
dark matter halo we must model its mass accretion history 
as well as the orbital evolution of the subsystems once they 
are accreted. For the first step, we rely on the the extended 
Press- Schechter (EPS) formalism to create merger histo- 
ries for each host system. We give a brief description of our 
EPS merger trees in §2.1. In §2.2 we discuss our model for 
the density structure of accreted halos and the host sys- 
tem and in §2.3 we describe our method for following the 
orbital evolution of each merged system. We show tests 
and examples of this model in §2.4. 

2.1. Merger Histories 



We track diffuse mass accretion and satellite halo acqui- 
sition of host systems by constructing merger histories us- 
ing the EPS method (Bond et al. 1991; Lacey & Cole 1993, 
LC93 hereafter). In particular, we employ the merger tree 
algorithm of Somcrville and Kolatt (1999, hereafter SK99). 
This allows us to generate a list of the masses and accre- 
tion redshifts of all subhalos greater than a given threshold 
mass that merged to form the host halo. We describe the 
method briefly here, and encourage the interested reader 
to consult LC93 and SK99 for further details. 

A merger tree that reproduces many of the results of N- 
body simulations can be constructed using only the linear 
power spectrum. For convenience, we express this in terms 
of a{M), the rms fluctuation amplitude on mass scale M 
at z = 0. As in LC93, let S{M) = (J^{M), AS = S{M) - 
S{M + AM), w{t) = 6c{t), and Sw = w{t) - w{t + At). 
Here 6c{t) is the linear overdensity for collapse at time 
t, associated with our choice of cosmology (see LC93 or 
White 1996). The probability that a halo of mass M, at 
time t, accreted an amount of mass associated with a step 
of AS', in a given time step implied by 5w is 



P{AS, 6w)d{AS) = 



dw 



27rA53/2 



exp 



-{6wf 
2AS 



d{AS). 



(1) 

Merger histories are constructed by starting at a chosen 
redshift and halo mass and stepping back in time with an 
appropriate time step. If the minimum mass of a progen- 
itor that we wish to track is Mmin, then SK99 tell us that 
each time step must be small in order to reproduce the con- 
ditional mass functions of EPS theory: 6w< ^Mrain{dS{M)/dM). 

We build merger trees by selecting progenitors at each 
time step according to equation (1) and treating events 
with AM < Mmin as diffuse mass accretion. At each step, 
we identify the most massive progenitor with the host halo 
and all less massive progenitors with accreted subhalos 
and we continue this process until the host mass falls be- 
low A'/,„in. In practice, we use a slightly modified version of 
the SK99 scheme. At each stage we demand that the num- 
ber of progenitor halos in the mass range we consider be 
close to the mean value. As discussed in BKW, this modi- 
fication considerably improves the agreement between the 
analytically predicted progenitor distribution and the nu- 
merically generated progenitor distribution. In what fol- 
lows we set Mmin = 10^ Mq. Our fiducial, z=0 host mass 
is 1.4 X lO'^^ M0, but we vary these choices in order test 
sensitivity to the host mass and redshift. 

2.2. Halo Density Structure 

Whether a merged system survives or is destroyed de- 
pends on the density structure of the subhalo and on the 
gravitational potential of the host system. Therefore, it 
is worthwhile to describe our assumptions about CDM 
density profiles in some detail. The size of a viriahzed 
dark matter halo can be quantified in terms of its virial 
mass A'/virj or cquivalently its virial radius i?vin or virial 
velocity V^^^ = GMvir/i?vir- The virial radius of a halo 
is deflned as the radius within which the mean density 
is equal to the virial overdensity Avir, multiplied by the 
mean matter density of the Universe pu, so that M^ir = 
4TTp\iAyir{z)RliJ3. The virial overdensity Avir, can be 
estimated using the spherical top-hat collapse approxima- 
tion and is generally a function of Om, ^a, and redshift 
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{e.g., Eke, Navarro, & Frcnk 1998). Wc compute Avir us- 
ing the fitting function of Bryan & Norman (1998). In the 
cosmology considered here, A^ij.{z = 0) ~ 337, and at high 
redshift Avir — > 178, approaching the standard cold dark 
matter {i.e., flu = 1) vahic. 

The gross structure of dark matter halos has been de- 
scribed by several analytic density profiles that have been 
proposed as good approximations to the results of high- 
resolution N-body simulations (Moore et al. 1999b; Power 
et al. 2003). In the interest of simplicity, wc choose 
to model all halos with the density profile proposed by 
Navarro, Frenk, & White (1997; hereafter NFW): 

^^"■^ " {r/n){l + r/ny 
For the NFW profile, the amount of mass contained within 
a radius r, is 

M(< r) = M™-7^ (3) 

where x = r/r^, g{y) = ln(l + y) — y/(l + y), and the 
concentration parameter is defined as Cyir = Rvir/rg. Re- 
stating equation (3) in terms of a circular velocity profile 
yields V^{r) = l/f;j.Cvir5'(a;)/.Tg(cvir). The maximum circu- 
lar velocity occurs at a radius rmax — 2.16rs, with a value 
F4,~0.216Klc™/ff(cvir). 

As a result of the study by Wechsler et al. (2002; W02 
hereafter) and several precursors {e.g., Zaroubi and Hoff- 
man 1993; NFW; Avila-Reese, Firmani, & Hernandez 1998; 
Bullock et al. 2001, hereafter BOl), we now understand 
that dark matter halo concentrations are determined al- 
most exclusively by their mass assembly histories. The 
gross picture advocated by W02 is that the rate at which 
a halo accretes mass determines how close to the center of 
the host halo the accreted mass is deposited. When the 
mass accretion rate is high, near equal mass mergers are 
very likely, and dynamical friction acts to deposit mass 
deep into the interior of the host. After an early period 
of rapid mass accretion, the central densities of halos re- 
main constant at a value proportional to the mean density 
of the Universe at the so-called "formation epoch" Zc, de- 
fined as the redshift when the relative mass accretion rate 
was similar to the rate of universal expansion (see W02 for 
details). For typical halos, this formation epoch occurred 
at a time when halos were roughly ~ 10 — 20% of their 
final masses. Additionally, W02 found that the scale ra- 
dius and central density of the best fit NFW profile remain 
practically constant after the initial phase of rapid accre- 
tion. After this time, the mass increase is slow, and as the 
virial radius of the halo grows, its concentration decreases 
as Cvir oc {1 + z)~^. 

These results lend support to BOl, who explained the 
observed trends with halo mass and redshift using a sim- 
ple, semi-analytic model that we adopt in this study. In 
the BOl model, halo concentrations Cvir(M, z), depend only 
on the value of cr(M) and the evolution of linear perturba- 
tions, 6{z)/S{z = 0). Specifically, the density of a halo of 
mass M is set by the density of the Universe at the time 
when systems of mass O.OIM were typically collapsing. 
The collapse epoch Zc, is defined by cr(O.OlM) = Sc{zc). 
Again, 6c{z) is the linear overdensity for collapse at red- 
shift z. Central densities determined in this manner con- 
nect well to the findings of W02. Halo density structure is 
set at a time of rapid accretion, when progenitor masses 



typically were Mprog ^ 0.1 M. Most of the mass in a halo 
at any given time is set by accretion events with subha- 
los of mass ^ 1/10 the host halo mass (c/. §4). Thus the 
period of rapid mass accretion involves objects of mass 
^ O.lA'/piog ^ O.OIM, and it is the collapse times and 
densities of these constituents that set the central density 
of the mass M halo. 

The BOl model reproduces N-body results for n = 1, 
n = 0.9, and power-law CDM models {e.g., Colin et al. 
2000b; BOl) as well the WDM simulations of Colin et 
al. (2000a) and Avila-Reese et al. (2001). However, we 
stress that N-body tests were restricted to the mass range 
~ 10^ — 10^** M0 because of the limited dynamic range 
of numerical experiments.'' Nevertheless, we use the BOl 
model to compute concentrations for halos with masses 
< 10^ M0. Our results for M< 10^ Mq may be regarded 
as a "best-guess" extrapolation of N-body results. 

Before proceeding, we mention an alternate prescription 
for assigning Cvir(M, z) proposed by Eke, Navarro, & Stein- 
metz (2001, ENS hereafter). ENS investigated the power 
spectrum dependence of the Cyir — Myir relation for sev- 
eral ACDM and WDM models. While the BOl and ENS 
recipes for Cvir(M, z) matched well for ACDM, the BOl 
model failed to reproduce the mass dependence seen in 
simulations by ENS for WDM halos with masses smaller 
than the "free-strc^aming" mass (see §3). The four WDM 
halos simulated by ENS with masses small enough to be 
appreciably affected by free-streaming all had Cvir values 
that were ~ 2a lower than the BOl model. Based on these 
data, ENS proposed a model in which halo collapse time 
depends not only on the amplitude of the power spec- 
trum a{M), but on an effective overdensity amplitude, 
(Tcff = -(j(M)dlncr(M)/dlnM. This results in a Cvir(M) 
relation that increases with mass for masses smaller than 
the truncation scale and decreases at larger masses as in 
ACDM. By defining an effective overdensity in this way, 
ENS were able to account for the low c^ir values observed 
in their WDM simulations and still reproduce the red- 
shift and mass dependence seen in ACDM simulations. 
The slope of the Cvir-Afvir of ENS is shallower than the 
slope predicted by the BOl relation, therefore the ENS 
model also leads to less concentrated halos at small mass 
{M^ 10^" Mq) even for identical input power spectra. 
This disparity grows larger when tilted and/or running 
spectra are considered, as in this paper. 

Unfortunately, the ENS model cannot be applied in the 
WDM cases we explore because in these models a{M) is 
very flat on scales smaller than the free-streaming mass 
and the ENS model breaks down when da{M) /dM be- 
comes very small. In the ENS model, WDM halos smaller 
than ^ 1% of the free-streaming mass never collapse be- 
cause (Tcs <C 1. In addition to this practical problem, the 
ENS predictions are not supported by the results of Avila- 
Reese et al. (2001) and Colin et al. ' (2000a). Using - 25 
halos, Avila-Reese et al. found WDM halo concentrations 
to be roughly constant with mass down to several orders of 
magnitude below the free-streaming scale, in accordance 
with the BOl model predictions. In light of these difficul- 
ties and the discordant results of different N-body studies, 

* Preliminary results from new simulation data show promising 
agreement with the BOl model all the way down to M ~ 10^ Mq 
(P. Colin, private communication). 
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we have not explored the implications of the ENS model 
in this work. This is not an indictment of the ENS model. 
Rather, the results of ENS highlight the uncertainty in 
assigning halo concentrations to low-mass systems, espe- 
cially with power spectra that vary rapidly with scale. Our 
choice of the BOl relation is a matter of pragmatism and 
represents a conservative choice in that halos are assigned 
the higher of the two predictions of Cvir at small mass. 
Lower Cvir values (in line with ENS expectations) would 
result in less substructure and larger deviations from the 
standard ACDM model than the predictions in §4. 

2.3. Orbital Evolution 

With the accretion history of the host halo in place, and 
with a recipe in hand that fixes the density structure of 
host and satellite halos, the next step is to track the orbital 
evolution of accreted systems. This is necessary in order 
to account for the effects of dynamical friction and mass 
loss due to tidal forces. These processes cause most of the 
accreted subhalos either to sink to the center of the host 
halo and become "centrally merged," or to lose most of 
their mass and be "tidally disrupted" and no longer iden- 
tifiable as distinct substructure. We model these effects 
using an improved version of the BKW technique, borrow- 
ing heavily from the dynamical evolution model proposed 
by Taylor and Babul (2001, TBOl hereafter; see also Tay- 
lor & Babul 2002) and the dynamical friction studies of 
Hashimoto, Funato, & Makino (2002, HFM02 hereafter) 
and Valenzuela & Klypin (2003; and Valenzuela & Klypin, 
in preparation). 

We denote the mass of an accreted subhalo as Mgat, its 
outer radius as -Rsat, and the accretion redshift as Zacc- We 
set the subhalo concentration to the median value given 
by the BOl model for this mass and redshift. Although 
initially set by the virial mass and radius of the in-falling 
subhalo, Msat and i?sat are allowed to evolve with time, 
as described in more detail below. We track the orbit of 
each subhalo in the potential of its host from the time of 
accretion iacc, until today [to ~ 13.5 Gyr in this cosmol- 
ogy) or until it is destroyed. The mass accretion history 
also yields the host halo mass at each time step. We fix 
the density profile of the host at each accretion time using 
the median BOl expectation for a halo of that mass. As 
we mentioned earlier, the scale radius and central density 
of the host remain approximately constant. 

For the purpose of tracking each subhalo orbit, we as- 
sume the host potential to be both spherically symmet- 
ric and static. We update the host halo profile using the 
BOl expectation at each accretion event, but hold it fixed 
while each orbit is integrated. While the approximation of 
a static host potential for each orbit is not ideal, it allows 
for an extremely simple prescription that significantly re- 
duces the computational expense of our study. Moreover, 
this approximation is not bereft of physical motivation. 
As we discussed above, halos observed in numerical simu- 
lations appear to form dense central regions early in their 
evolution after which their scale radii and central densities 
remain roughly fixed with time.^ Additionally, we have 
run test examples that include an evolving halo potential 

^ The exception to this is the case of a late-time merger of halos of 
comparable mass in which case the central densities and scale radii 
of the participating halos may change considerably (W02). 
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Fig. 2. — Input orbital circularity distribution of initially in-falling 
substructure (dashed) shown along with the circularity distribution 
of the final surviving population of (n = 1) LCDM subhalos at z = 
(solid). For reference, the thin dotted line shows the circularity 
distribution of surviving substructure measured by Ghigna et al. 
(1998) in their N-body simulations. 



(set by the results of W02) and find that this addition has 
a negligible effect on the statistical properties of substruc- 
ture that we are concerned with here. 

Upon accretion onto the host, each subhalo is assigned 
an initial orbital energy based on the range of binding 
energies observed in numerical simulations (K99; A. V. 
Kravtsov 2002, private communication). We place each 
satellite halo on an initial orbit of energy equal to the en- 
ergy of a circular orbit of radius i?ciic = ^^vir, where i?vir 
is the virial radius of the host at the time of accretion 
and rj is drawn randomly from a uniform distribution on 
the interval [0.4,0.75]. We assign each satellite an initial 
specific angular momentum J ~ eJcirc, where Jdrc is the 
specific angular momentum of the aforementioned circu- 
lar orbit and e is known as the "orbital circularity." Past 
studies drew e from a uniform distribution on the inter- 
val [0.1, 1] (BKW) to match the circularity distribution of 
surviving subhalos in simulations reported by Ghigna et 
al. (1998). However, the orbits of surviving halos are bi- 
ased relative to the orbits of all accreted systems because 
subhalos on radial orbits are preferentially destroyed. We 
find that we better match the Ghigna et al. (1998) result 
for surviving satellites if we draw the initial e from the 
simple, piecewise-linear distribution depicted in Figure 2. 
The initial radial position of each satellite halo is set to 
^init = ^circ and for all non-circular orbits, we set the 
subhalo to be initially in- falling so that dR/dt < 0. 

To calculate the trajectories of subhalos, we treat them 
as point masses under the infiuence of the NEW gravita- 
tional potential of the host halo. We model orbital decay 
by dynamical friction using the Chandrasekhar formula 
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(Chandrasekhar 1943). The Chandrasekhar formula was 
derived in the context of a highly idealized situation; how- 
ever, numerical studies indicate that this approximate re- 
lation can be applied more generally (Valenzuela & Klypin 
2003 have performed a new test that supports the use of 
this approximation). Using the Chandrasekhar approxi- 
mation, there is a frictional force exerted on the subhalo 
that points opposite to the subhalo velocity: 



^DF 



47rln(A)G2M,i,p(r) 



erf(X) 



2X 



: exp(-X2) 



(4) 

In equation (4), ln(A) is the Coulomb logarithm, r is the 
radial position of the orbiting satellite, and p{r) is the 
density of the host halo at the satellite radius. The quan- 
tity Vorh is the orbital speed of the satellite halo and 
X = K)rb/V2(7^, where a is the one-dimcnsional veloc- 
ity dispersion of particles in the host halo. For an NFW 
profile, the one- dimensional velocity dispersion can be de- 
termined using the Jeans equation. Assuming an isotropic 
velocity dispersion tensor, 

a {x = r/rs) = Kir"? + x) 







-dx' . 



(5) 

We find the following approximation useful and accurate 
to 1% for x = 0.01 - 100: 

1 4393a;°'^^^ 

" ^-TTH756^- 

There has been much debate on the appropriate way 
to assign the Coulomb logarithm in Eq. (4). Dynami- 
cal friction is caused by the scattering of background par- 
ticles into an overdense "wake" that trails the orbiting 
body and tugs back on the scatterer. The Coulomb log- 
arithm is interpreted as ln(6max/^min)5 where 6max is the 
maximum relevant impact parameter at which background 
particles are scattered into the wake and 5niin is the min- 
imum relevant impact parameter. A common approxima- 
tion is to choose a constant value of the Coulomb loga- 
rithm (perhaps by calibrating to the results of numerical 
experiments as in TBOl), but some studies indicate that 
this approach significantly underestimates the dynamical 
friction timescale when tested against N-body simulations 
{e.g., Colpi, Mayer, & Governato 1999; HFM02). Moti- 
vated by the results of HFM02 and Valenzuela & Klypin 
(in preparation), we allow the Coulomb logarithm to evolve 
with time and set 6max = ^(0) where r{t) is the radial po- 
sition of the orbiting subhalo. We assign the minimum 
impact parameter according to the prescription of White 
(1976) and integrate the effect of encounters with back- 
ground particles over the density profile of the subhalo. 
Repeating this calculation for an NFW halo yields 



ln(A) = In 



1 



at 



where 

/(X,at) = Xl 

JO Jxb 



I{Xsim), 



(7) 



X^ 



--dx 



dxb, (8) 



Xsa.t = -Rsat/^'rS ^'^d ''s'^* is the NFW scale radius of the 
satellite. The integral I{y) is well-approximated by the 



following function, which is accurate to 1% for 0.1 < y < 
100: 

0.10947?/-"«" 



[1 + 0.90055?/ 



1.099 . 



- 0.03568J/ 



1.189 . 



0.064032/1 



(9) 

As the satellite orbits within the host potential, it is 
stripped of mass by the tidal forces that it experiences. 
First, we estimate the instantaneous tidal radius of the 
subhalo rt, at each point along its orbit. In the limit that 
the satellite is much smaller than the host, the tidal radius 
is given by the solution to the equation (von Hoerner 1957; 
King 1962) 

3 M,at(< rt)/Mhost(< r) 3 



't — 



2 + a;2ij3/GMhost(< r) - ainMhost(< r)/d\n, 



(10) 

where r is the radial position of the satellite, Mhost(< r) is 
the host's mass contained within this radius [cf. Eq. (3)], 
.^sat(< ''t) is the satellite's mass contained within rt, and 
oj is the instantaneous angular speed of the satellite. Equa- 
tion (10) is merely an estimate of the satellite's tidal limit. 
For a satellite on a circular orbit, it represents the distance 
from the satellite center to the point along the line con- 
necting the satellite and the host halo center where the 
tidal force on a test particle just balances the attractive 
force of the satellite. In reality, the tidal limit of a satellite 
cannot be represented by a spherical surface: some parti- 
cles within rt will be unbound while others without rt may 
be bound. Nevertheless, TBOl showed that this can serve 
as a very useful approximation. 

As the tidal radius shrinks, unbound mass in the periph- 
ery is stripped. Tidal forces are strongest, and rt smallest, 
when the orbit reaches pericenter; however, all of the mass 
outside of rt is not stripped instantaneously at each peri- 
center passage. Rather, mass is gradually lost from the 
satellite on a timescale set by the orbital energy of the lib- 
erated particles. Johnston (1998), found that the typical 
energy scale of tidally stripped debris is set by the change 
in the host halo potential on the length scale of the or- 
biting satellite, e w rtd$host(?')/dr. Particles on circular 
orbits of energy E and E±e move a distance rt relative to 
each other on a timescale of order the orbital period, T. As 
such, we may expect T to be the relevant timescale for tidal 
mass stripping. TBOl used this timescale in their model 
to reproduce the results of several idealized N-body exper- 
iments. Following TBOl, we model satellite mass loss by 
dividing the orbit into discrete time steps of size dt <C T. 
At each step, we remove an amount of mass 

Sm = M,at(> rt)j, (11) 

where Msat(> rt) is the satellite's mass exterior to rt. 

As a subhalo loses mass due to tidal stripping, we as- 
sume that its density profile is unmodified within its outer 
radius -Rsat- Rather than identify i?sat with the tidal ra- 
dius (which does not vary monotonically with time), we 
set its value by determining the radius within which the 
mass profile retains the appropriate bound mass using Eq. 
(3). We fix the scale radius of the subhalo rg**', at the 
value defined at the epoch of accretion. 

Our approximations for dynamical friction and tidal strip- 
ping are least accurate when the mass of the satellite is not 
very small compared to the mass of the host. However, as 
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Fdf oc MjHj^., it is in precisely these cases that we expect 
the satellite to merge quickly with the host and no longer 
be identifiable as distinct substructure. As such, the pre- 
cise dynamics should not have a significant effect upon 
our results in these cases, particularly because our main 
predictions involve low-mass substructure. However, more 
detailed modeling will be important for investigations that 
focus on more massive substructures, for example, explo- 
rations that use disk thickening as a test of the ACDM 
cosmological model {e.g., Font et al. 2001). 

The final ingredients for our semi-analytic model of halo 
substructure are the criteria for declaring subhalos to be 
tidally disrupted and centrally merged. Let r^^x — ^.IQrl'^^ 
be the radius at which the subhalo's initial velocity profile 
attains its maximum, and Msat(< ^max) be the mass of 
the satellite originally contained within the radius r^ax- 
We declare a subhalo to be centrally merged with the host 
if its radial position relative to the center of the host be- 
comes smaller than r^*^. We declare a satellite tidally 
disrupted if the mass of the satellite becomes less than 
Msnt{< '"max)- This Criterion is partially motivated by the 
numerical study of H03, who find that NFW subhalos are 
completely tidally destroyed shortly after rt becomes less 
than r^*x- Of course the distinction between centrally 
merged and tidally destroyed satellites is somewhat arbi- 
trary as subhalos are typically severely tidally disrupted 
as they approach the center of the host potential. Fortu- 
nately, for the issues we explore here, the precise nature of 
a satellite halo's destruction is not important. We discuss 
this issue further in a forthcoming extension of our work 
(A. Zentner & J. Bullock, in preparation). 

In reporting results concerning the velocity function of 
substructure, we invoke a further modification. H03 noted 
that subhalos that experienced significant tidal stripping 




Fig. 3. — Orbital evolution for three sets of subhalo parameters: 
ACt = 10» Mo , Cvir = 15 (solid); M^^^ = 10* Mq , c^i, = 7.5 
(dashed); and M°^^ = 5 x 10^ Mq , Cvir = 15 (short-dashed). Initial 
orbital parameters and host mass properties are fixed, as described 
in the text. The top panel shows the radial evolution in units of 
the initial radius as a function of time. The bottom panel shows 
the mass of each system as a function of time. Lines that terminate 
represent subhalo destruction at the end point (see text). 



satellite was started on the same initial orbit, e = 77 = 0.5, 



suffered not only mass loss at radu > rt, but mass redis- but the satellite properties ^were varied: M^^^ = 10^ Mq, 
tribution in their central regions, at radii smaller than rt. 
To account for this, we determine whether or not the tidal 
radius of each surviving subhalo was ever less than r^^x- 
If so, we follow the prescription of II03 to account for mass 
redistribution and scale the maximum circular velocity of 
the satellite via 



-^final 



Affinal \ 



ini 
max 



tial 



(12) 



where X^ax is the maximum circular velocity of the satel- 
lite according to its initial density profile, M^^"^' is its final 
mass, and Mg™'"*' is its initial mass before being tidally 
stripped. In practice, this rescaling has a fairly small ef- 
fect on our velocity functions. Roughly ~ 30% of surviving 
halos meet this condition for rt. For those halos that do 
experience this kind of mass loss, the typical reduction in 

Knax is < 25%. 

We are currently in the process of checking this model 
against idealized N-body experiments designed to mimic 
the type of orbital histories that we encounter here (J. 
Bullock, K. Johnston, & A. Zentner, in preparation). Pre- 
liminary results show promising agreement. 

2.4. Tests and Examples 

Figure 3 shows three example calculations of subhalo 
trajectories aimed at demonstrating how various factors 
affect the orbital evolution of a satellite system. Each 



Cvii- = 15 (solid); M^^t = 10^ Mq, Cvir = 7.5 (dashed); and 
M°^^ = 5 X 10^ Mq, Cvir = 15 (short-dashed). The up- 
per and lower panels depict the evolution of orbital radius 
and mass of the subhalo respectively. The accretion time 
was set at 8 Gyr in the past, a = {1 + z)~^ ~ 0.45 for 
this cosmology. The host halo parameters were chosen to 
match reasonable expectations for a Milky Way-sized pro- 
genitor at that time: Mhost = 5x 10^^ Mq (i?vir — llOkpc) 
and Cvir = 6. While the subhalo represented by the solid 
line experiences gradual tidal mass loss and slight orbital 
decay as a result of dynamical friction, its core survives 
for the full time period. The less concentrated subhalo 
(dashed) is more strongly affected by tides, and is com- 
pletely disrupted ~ 3.5 Gyr after being incorporated into 
the host. (Although not shown, a similar effect is seen 
if the host halo concentration is increased and the sub- 
halo concentration is held fixed.) In the case of the mas- 
sive subhalo, dynamical friction causes the orbit to decay 
more quickly and the subhalo experiences more frequent 
pericenter passages. Consequently, disruption occurs ~ 6 
Gyr after accretion. Notice that because the stripping 
process is gradual (unless orbits are very radial) and the 
timescales involved are of order ~ Gyr, the accretion time 
is also important in determining survival probability. If 
any of these subhalos were were accreted more recently, 
their chance of survival to the present day would increase 
accordingly. The combination of factors illustrated here 
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Fig. 4. — The velocity functions of progenitor and surviving sub- 
halo populations derived using our fiducial (n = 1, erg = 0.95) 
ACDM cosmology and a 200-halo ensemble of 1.4 X 10^'^ Mq systems 
at 2 = 0. Shown are all accreted halos (dashed), and the fraction of 
those that are tidally destroyed (short-dashed) and centrally merged 
(dotted). The solid line shows the surviving population of subhalos 
at z = and, for comparison, the thin dashed line shows the sur- 
viving population derived by K99 using N-body simulations. The 
error bars represent the sample variance. 



— accretion time, satellite mass, and the relative concen- 
trations of host and satellite — will be important in later 
sections for understanding the factors that set the subhalo 
population from one cosmology to the next. 

Figure 4 shows the ensemble-averaged, cumulative ve- 
locity function for the substructure population of Milky 
Way-like host halos computed in our standard ACDM cos- 
mology. The host properties a,t z = are M^h- — 1-4 x 10^^ 
Mq, Cvir — 13.9, and Knax — 187 km s~^. The lines rep- 
resent the means of 200 merger tree realizations, and the 
error bars represent the sample variances over these re- 
alizations. In particular, the thick solid line shows the 
surviving subhalo population at z = 0. For comparison, 
the thin dashed line is the best-fit velocity function re- 
ported by K99 based on an analysis of substructure in 
ACDM halos. This line is plotted over the range that 
their resolution and sample size allowed them to probe. 
The apparent agreement between our semi-analytic model 
and the N-body result is excellent, and lends confidence in 
our ability to apply this model to different power spectra. 

The radial distribution of substructure at z = for the 
same ensemble of halos is shown in Figure 5. Open cir- 
cles show the differential number density profile of subha- 
los with Afsat > 10^ Mq normalized relative to the total, 
volume-averaged number density of subhalos within i?vir 
that meet the same mass requirement. The solid pen- 
tagons show the same quantity for more massive subha- 
los, Msat > 6 X 10^ Mq. The line shows the NFW dark 
matter profile for the host system normalized relative to 




r (kpc) 



Fig. 5. — The radial number density profile of substructure derived 
from 200 model realizations of a Mvir = 1.4 X lO^'^ Mq host halo 
at 2 = in our fiducial (n = 1, erg = 0.95) ACDM cosmology. The 
open circles show the number density of subhalos with A/sat > 10® 
Mq divided by the average number density of systems meeting this 
mass threshold within the virial radius of the host system. The 
points reflect the radial profile averaged over all realizations, and 
the error bars reflect the sample variance. Solid pentagons show the 
same result for A/sat > 6 X 10® Mq subhalos. The variance (not 
shown) is significantly larger for the higher mass threshold because 
there are significantly fewer such systems in each host. For reference, 
the solid line shows the NFW density profile of the host at 2 = 0. 
The virial radius for a host halo of this size is Hvir — 285 kpc 
and the typical NFW scale radius is Ts ~ 20 kpc. We do not plot 
predictions beyond r = 0.75-Rvir — 215 kpc because this is the 
maximum circular radius we assign to in-falling, bound systems. 



the average (virial) density within the halo. Observe that 
the subhalo profile traces the dark matter profile at large 
radius (p oc r~^), but fiattens towards the center as a con- 
sequence of tidal disruption. This result agrees well with 
that presented in Figure 3 of Colin et al. (1999). Us- 
ing an N-body analysis of a cluster-sized host, Colin et al. 
(1999) showed that the number density of systems with 
Msat greater than 0.04% of the host mass traces the back- 
ground halo profile at large radius, begins to flatten at 
r ~ 0.2i?vir, and is roughly a factor of 5 below the back- 
ground at r ~ 0.07i?vir (their innermost point). ^ The solid 
pentagons in Figure 5 correspond to the same mass frac- 
tion relative to the host. Notice that at r = 0.07i?vir — 20 
kpc, the factor of '--^ 5 mismatch is reproduced. Ghigna 
et al. (1998) observed the same qualitative behavior for 
subhalos in a standard CDM simulation of a cluster-size 
halo. Chen et al. (2003) have measured the substruc- 
ture profile using a high-resolution galaxy-size halo with 
Msat^ 0.0015%Mhost, corresponding to subhalos interme- 
diate in mass between those represented by the open cir- 

^ We quote results relative to iJvir and Afjjost because the host halo 
in Colin et al. (1999) is significantly more massive than the halos 
that we consider. 
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cles and solid pentagons in Figure 5. Chen et al. (2003) 
similarly find core behavior setting in at a radius of ^ 30 
kpc, but also find a stronger overall suppression in sub- 
structure counts within 70 kpc. Our results suggest 
that some of the observed suppression may be caused by 
overmerging in the central regions of their simulated halos. 
Ongoing studies by other workers lead to similar conclu- 
sions (J. Taylor, private commimication) . Only the next 
generation of numerical simulations can reliably test this. 
That we produce a reasonable approximation to the num- 
ber density profile of substructure is an indication of the 
soundness of our model. 

3. MODEL POWER SPECTRA 

The initial power spectrum of density fluctuations is 
conventionally written as an approximate power law in 
wavenumber k, P{k) oc k", corresponding to a variance per 
logarithmic interval in wavenumber of A'^{k) = fc^P(fc)/27r^. 
If the fluctuations were seeded during an early inflation- 
ary stage, as is commonly supposed, then the initial spec- 
trum is likely to be nearly scale- invariant, with n ~ 1. 
Any deviation from power law behavior, or "running" of 
the power law index with scale is likewise expected to be 
smaU, |dn/dlnfc| < 0.01 (Kosowsky & Turner 1995). In 
addition to these theoretical prejudices, large-scale obser- 
vations of galaxy clustering and CMB anisotropy seem to 
favor nearly scale-invariant models that can be parameter- 
ized in this way. In this paper we explore the effects on 
halo substructure of taking n ^ 1 and allowing for scale- 
dependence in the power law index and more dramatic 
features in the power spectrum. In this section we give a 
brief description of the power spectra that we explore. 

Table 1 summarizes the relevant features of our example 
power spectra. The second and third columns list the pri- 
mordial spectral index evaluated at the pivot scale of the 
COBE measurements fccoBE ~ 0.0023 h Mpc~^, and the 
running of the spectral index. ^ We neglect any variation in 
the running with scale and explicitly set d^n(fc) /d(ln fc)^ = 
0. Except for the running index (RI) case, we normal- 
ize all models to the COBE measurements of the CMB 
anisotropy using the fitting formulae of Bunn, Liddle, & 
White (1996; also Bimn & White 1997). The fourth col- 
umn of Table 1 gives the implied value of as ■ We calculate 
spectra using the transfer functions of Eisenstein & Hu 
(1999). In Figure 6 we illustrate the implied a{M) for 
these models. 

Many of the spectra listed in Table 1 are motivated by 
particular models of inflation. We invoke an inverse power 
law potential that gives rise to a mild tilt n ~ 0.94, as 
well as a model in which the inflaton has a logarithmically 
running mass and which can give rise to significant tilt and 
running for natural parameter choices (Stewart 1997a, b; 
Covi & Lyth 1999; Covi, Lyth, & Roszkowski 1999; Covi, 
Lyth, & Melchiorri 2003). We employ specific inflationary 
potentials mainly as a conceptual follow up to ZB02, which 
highlighted the fact that various levels of tilt may occur 
naturally within the paradigm of inflation and that n = 1 
is not demanded by this paradigm. For the purposes of 
this paper, one may regard our choices simply as spanning 

Wc use the definition of running employed by Spcrgcl et al. (2003) 
rather than that given in, for instance, Kosowsky & Turner (1995). 
These definitions diflfer by a factor of two. 



a range of observationally viable input power spectra. The 
values of tilt and (jg that we consider range from n ~ 0.84 
with ag = 0.65 to n = 1 and erg = 0.95. The model with 
fTg = 0.75 was specifically chosen to match galaxy central 
densities, as described in ZB02. We also explore the best- 
fit, running-index model of the WMAP team (Spergel et 
al. 2003), with dn/dlnfc = -0.03. We refer to this as the 
"running index model" or "RI model." Note that Spergel 
et al. (2003) quote a value of n = 0.93 evaluated at k = 
0.05 Mpc^^. The value listed in Table 1 is larger because 
we quote it at a smaller wavenumber, k = fccOBE- 

In addition to tilted ACDM models, we consider spec- 
tra with abrupt reductions in power on small scales. In 
the "broken scale-invariancc" (BSI) example, we adopt an 
idealized inflation model introduced by Starobinsky (1992) 
that exhibits the most rapid drop in power possible for a 
single field model. Kamionkowski & Liddle (2000) studied 
this type of model as a way to mitigate the dwarf satellite 
problem, but our choice of parameters is slightly different 
from theirs (see ZB02). 

We also consider WDM scenarios in which the primor- 
dial power spectra are scale-invariant but small-scale fluc- 
tuations are subsequently filtered by free-streaming. The 
free-streaming scale is set by the primordial velocity dis- 
persion of the warm particles. In the canonical case of 
a "neutrino-like," thermal relic with two internal degrees 
of freedom, the free-streaming scale can be expressed in 
terms of the warm particle mass mw and relic abundance, 



i?/ ~0.11 





1/3 




0.15 




IkevJ 



-4/3 



Mpc. (13) 



We calculate WDM spectra assuming the same flat cosmol- 
ogy with f2M = ^^w + = 0.3, and use the approximate 
WDM transfer function given by Bardeen et al. (1986), 
P{k) = exp[-kRf - (kR.jf]PcBM{k). 

Several studies have placed approximate constraints on 
WDM masses based on either the argument that there 
must be enough power on small scales to reionize the Uni- 
verse at sufficiently high redshift (zj-cJ^ 6) or by probing 
the power spectrum on small scales directly with the Ly- 
a forest (Barkana, Haiman, & Ostriker 2001; Narayan et 
al. 2000). These authors essentially find that mw^ 0.75 
keV assuming a neutrino-like thermal relic; however, this 
constraint may be significantly more restrictive if measure- 
ments of Zre ~ 17 by the WMAP collaboration (Kogut et 
al. 2003; Spergel et al. 2003) are confirmed (Somerville, 
Bullock, & Livio 2003). As such, we consider three illus- 
trative examples in what follows, mw = 0.75 keV, 1.5 keV, 
and 3.0 keV. The corresponding "free-streaming" masses, 
below which the fiuctuation amplitudes are suppressed, are 
listed in Table 1. 

4. RESULTS 

4.1. Accretion Histories 

Our first results concern the merger histories of halos 
that are approximately Milky Way-sized, with M^ir = 
1.4 X 10^2 jvIq at z = 0. For the n = 1, ACDM model, 
wc present results based on 200 realizations. For all other 
models, our findings are based on 50 model realizations. 
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Fig. 6. — The 2 = rms overdensity as a function of mass scale for several of our adopted power spectra (see §3 and Table 1 for more 
details). In the left panel, we exhibit spectra that deviate from the standard n = 1, scale-invariant model. The models shown in this panel are 
standard n = 1 (solid), a broken scale-invariant model (dotted), eg = 0.84 & n ~ 0.94 (short-dashed), trg = 0.65 & n ~ 0.84 (long-dashed), 
erg = 0.75 & n ~ 0.90 (dot-long-dashed), and a model based on the results of the WMAP team with n ~ 1.03, dn/dlnfc = —0.03, and 
erg = 0.84 (dot-short-dashed). The inflation models that inspire these examples are described in ZB02. In the right panel, we show several 
warm dark matter power spectra. More precisely, we depict spectra implied by warm particle masses of mw = 3.0 keV (short-dashed), 
mw = 1-5 keV (long-dashed), and rrtw = 0.75 keV (dotted) along side the standard n = 1, ACDM spectrum (solid). 



Table 1 



Summary of power spectra properties 



Model Description 


"(fccOBE) 


dn(fc)/d Infc 


0-8 


comments 


Scale-invariant 


1.00 


0.000 


0.95 




Inverted power law inflation 


0.94 


-0.002 


0.83 




Running-mass inflation I 


0.84 


-0.008 


0.65 


see Stewart 1997a, b 


Running-mass inflation II 


0.90 


-0.002 


0.75 




WMAP best-fit running index (RI) model 


1.03 


-0.03 


0.84 


WMAP best fit, see Spergel et al. 2003 


Broken scale-invariant inflation (BSI) 


1.00 


0.000 


0.97 


exhibits sharp decline in power at 1 h Mpc~^, 










power suppressed for 10^" Mq 


Warm Dark Matter, mw = 3.0 keV 


1.00 


0.000 


0.95 


Mf ~ 8.3 X 10* M0 


Warm Dark Matter, mw = 1-5 keV 


1.00 


0.000 


0.95 


Mf ~ 1.3 X 10^0 M0 


Warm Dark Matter, mw = 0.75 keV 


1.00 


0.000 


0.94 


Mf ~ 2.1 X 10" Mq 



Note. — Column (1) gives a brief description of the inflation or warm dark matter model used to predict the power spectrum. In the text, 
we distinguish the first five models by their tilts and/or their values of erg. We label the warm dark matter models by the warm particle mass. 
Columns (2) and (3) give the tilt n(fccoBE) on the pivot scale of the COBE data fccoBE ~ 0.0023 h Mpc~^, and the running of the spectral 
index dn(fc)/dlnA;, respectively. We have explicitly assumed the "running-of-running" to be small and taken d^n(fe)/d(lnA;)^ = 0. Column (4) 
contains the values of erg implied by the tilt or warm particle mass, the COBE normalization, and our fiducial cosmological parameters except 
in the case of the WMAP best-fit running index (RI) model, in which case the value of erg reflects their best-fit normalization. 
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Fig. 7. — Fraction of final host mass accreted in subhalos of mass 
Msat as a function of Msat- The final host mass is 1.4 X 10^^ Mq. 
The results for several input power spectra are shown. 



Figure 7 focuses on the the mass distribution of accreted 
halos, integrated over the entire merger history of the host. 
We plot d//dlog(Afsat), the fraction of mass in the final 
halo that was accreted in subhalos of a given mass per 
logarithmic interval in subhalo mass, Mgat- Observe that 
the mass fraction accreted in subhalos of a given mass is 
relatively insensitive to the shape of the power spectrum. 
Although similarity from model to model may be some- 
what surprising at first, it follows directly from repeated 
application of Equation (1). In particular, the shape of 
the progenitor distribution for Mgat ^ -/Vfhost must follow 

1/2 

d//dlog(Msat) oc Mgj^j , and the turnover occurs because 
mass conservation suppresses the number of major merg- 
ers. The shape shown in Figure 7 and its insensitivity to 
the power spectrum is discussed in detail by LC93. 

While the total mass function of accreted substructure is 
relatively independent of the power spectrum, the merger 
histories themselves are not. In models with less power on 
galaxy scales, halos assemble their mass later and experi- 
ence more recent mergers and disruption events. We show 
an example of this in Figure 8. Here we plot the average 
accretion rate of subhalos with Msat > 10* Mq for host 
halos in the standard n — 1, ACDM model, the RI model, 
and our lowest normalization case {n = 0.84, as = 0.65). 
The total accretion rate is divided in two pieces: dashed 
lines show those subhalos that are eventually destroyed 
and solid lines show the accretion times of subhalos that 
survive until z = 0. For the standard (n = 1, cg = 0.95) 
case, the event rate peaks sharply about ~ 12 Gyr in the 
past, while the erg = 0.65 case has a broader distribution, 
peaking later at ^ 9 Gyr ago, and with a long tail of ac- 
cretion events extending towards the present day. 

The shift in accretion times in models with less small- 




5 10 
(lookback) time of accretion (Gyr) 

Fig. 8. — Accretion rate averages, dN/dt (Gyr~^), for merged 
halos more massive than 10* and host halos of mass 1.4 X lO^'^ 
Mq at 2 = 0. Three different power spectra are shown: n = 1 
(upper); RI model (middle); and n = 0.84, erg = 0.65 (bottom). 
Dashed lines show objects that are destined to be destroyed, either 
by tidal disruption or central merging, and solid lines show subhalos 
that survive until 2 = 0. The n = 1 results are based on 200 EPS 
realizations while the others are based on 50 realizations. 



scale power plays an important role in regulating the num- 
ber of surviving subhalos. As we discussed in relation to 
Figure 3, a finite amount of time is required for an orbit 
to decay or for a system to become unbound and in many 
cases the longer a subhalo orbits in the background poten- 
tial, the more probable its disruption becomes. The later 
accretion times in models with less power partially com- 
pensate for the fact that subhalos in these models are less 
centrally concentrated and more susceptible to disruption 
at each pericenter passage. Particular results for substruc- 
ture populations in each model are given in the following 
subsections. 

That we expect a characteristic merger/disruption phase 
in each halo's past is intriguing, as this phase is approxi- 
mately coincident with the estimated ages of galactic thick 
disks, itrf ~ 8 - 10 Gyr {e.g., Quillen & Garnett 2000 for 
the Milky Way), which seem to be ubiquitous and roughly 
coeval (Dalcanton & Bernstein 2002). In this context, the 
age distributions of thick disks might serve as a test of this 
characteristic accretion time, which varies as a function of 
normalization and cosmology. We reiterate that the look- 
back times shown for the dashed lines in Figure 8 are the 
times that the subhalos were accreted. The distributions 
of central merger rates and tidal destruction rates peak at 
slightly more recent times and their widths are broader, 
with longer tails towards the present epoch. 

It is interesting to note that the surviving halos in Figure 
8 represent a distinctly different population of objects than 
the destroyed systems — they tend to have been accreted 
more recently. We are inclined to speculate that the star 
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Fig. 9. — The cumulative mass function A'^(> M), of surviving 
subhalos computed for an ensemble of host halos of mass A/jjost = 
1.4 X 10^'^ Mq at z = 0. The lines show the means computed over 
all realizations. The different line types relate to different models 
following the convention in the left panel of Figure 6. The n = 
1 results are derived from 200 realizations and the results for the 
other models are based on 50 realizations. The error bars show the 
variance over the n = 1 realizations (upper) and BSI realizations 
(lower). 



formation histories of galaxies that were destroyed after 
being accreted could be distinctly different from those of 
the surviving (dwarf satellite) galaxies as well. This may 
have implications for understanding whether the stellar 
halo of our Galaxy formed from disrupted dwarfs or some 
other process. While the global structure of the stellar 
halo seems consistent with the disruption theory (Bullock, 
Kravtsov, & Weinberg 2001), the element ratios of stellar 
halo stars and stars in dwarf galaxies are not consistent 
with a common history of chemical evolution (Shetrone, 
Cote, & Sargent 2001). The results shown in Figure 8 
provide general motivation to model dwarf galaxy evolu- 
tion and Milky Way formation in a cosmological context. 

4.2. Mass and Velocity Functions 

We present our results on CDM halo substructure be- 
ginning with the abundance of satellites in Milky Way- 
like galaxies. We plot the mass function of subhalos iV(> 
Msat), or the number of subhalos with mass greater than 
Msat as a function of Mgat , for each of our models in Figure 
9. The host halo mass is again fixed at 1.4 x 10^^ Mq at 
z = 0. From this figure, we see that even in the signif- 
icantly tilted, low- normalization model (erg — 0.65), the 
number of satellite halos with mass greater than 10^ Mq 
is roughly equal to that in the standard n = 1 model. The 
systematic differences between models are small compared 
to the scatter. The suppression is weak because several 
competing effects tend to compensate for the reduced con- 
centrations of the subhalos. In tilted models with reduced 
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Fig. 10. — The cumulative velocity function of subhalos in a host 
of fixed mass, M^ost = 1.4 X lO-''^ Mq. The lines represent averages 
over 200 merger histories for n = 1 and 50 merger histories for the 
other models. The error bars represent the dispersion amongst these 
realizations. The models are n = 1 (solid), eg ~ 0.83 (dashed), RI 
model (dot-short-dashed), erg ~ 0.65 (dot-long-dashed), and BSI 
(dotted). 



small-scale power, subhalos are typically accreted at later 
times. In addition, host halos are less concentrated and 
correspondlingly less capable of disrupting their satellites. 

In contrast, the BSI model shows a substantial decrease 
(a factor of ^ 3) in the number of surviving satellite halos 
at fixed mass. The reason for the dramatic reduction in 
this case is easy to understand. First, power is reduced 
only on scales smaller than a critical scale around ~ 10^" 
M0 (c/. Figure 6) and so, the concentration and accretion 
history of the ~ 10^^ Mq host halo are minimally altered 
while the concentrations of the small subhalos are drasti- 
cally reduced (see ZB02). In other words, the host halo has 
a density structure similar the n — 1 model host and is just 
as capable of tidally disrupting satellites, but the satellites 
are significantly more susceptible to disruption. A second 
difference is that galaxy-size halos in the BSI model, unlike 
the tilted models, accrete ~ 40% fewer low-mass 10'' 
Mq ) halos over their lifetimes, and this further widens 
the disparity between the BSI and tilted-ACDM models. 

It is conventional to discuss the substructure popula- 
tion of Milky Way-like halos in terms of the velocity func- 
tion. In Figure 10, we show our results for the cumula- 
tive velocity functions of subhalos for a fixed host mass 
of Afhost = 1.4 X 10^^ Mq. Notice that the velocity func- 
tions show a stronger trend with power spectrum than 
the mass functions (Figure 9), but the effect is still mod- 
est compared to the statistical scatter. For the most ex- 
treme tilted model, the total number of subhalos with 
Knax^ 10 km s~^ is only a factor of ~ 2 lower than in 
the standard, scale-invariant case. In the case of the tilted 



14 



Andrew R. Zentner & James S. Bullock 



100 



10 



A 




n=l, CTa = 0.9b 
c7g-0.83 

d:i cliik 0.03 Rl 
(7„=0.65 



10 



20 

V,, (km s-i) 



40 



60 



Fig. 11. — The cumulative velocity function of subhalos in a 
host of fixed maximum velocity Knax — 187 km s~^. The fines 
represent averages over 200 merger history reafizations for n = 1 
and 50 realizations for all other models. The error bars represent 
the dispersion in these realizations. The different models are as in 
Figure 10. 



models, the reduction in the velocity function is largely 
due to the fact that the subhalos are less concentrated, so 
the I^nax values are correspondingly smaller for fixed halo 
masses [cf. Eq. (3) and the discussion that follows]. 

This effect is illustrated explicitly in Figure 11 where, 
rather than fixing the host mass at z = 0, we have fixed 
its maximum circular velocity at T^nax = 187 km s~^, the 
value of a typical n = 1, A/host = 1-4 x 10^^ Mq halo at 
z = 0. Normalizing our host halos by l^nax rather than 
mass is perhaps a more reasonable choice because Knax 
is more closely related to observations.^ Models with less 
galactic-scale power require a more massive host in or- 
der to obtain the same value of T^nax, and their veloc- 
ity functions shift correspondingly. For example, a host 
with Vniax — 187 km s^^ in the erg = 0.65 model requires 
Mhost ~ 2.2 X 10^^ M0. With this adjustment, the ve- 
locity functions of the various tilted models are now very 
similar. Again, the BSI case is different from the tilted 
models because the relative shift in the V^ax--^vir relation 
changes abruptly with mass scale. It is also encouraging 
that our model BSI velocity function agrees well with the 
N-body results of Colin et al. (2000) for a similar type of 
truncated power spectrum (see their Rf — 0.1 Mpc model. 
Fig. 2). 

Another convenient way to quantify the substructure 
content of halos is the fraction of mass in subhalos less 

* Vmax = 187 km s~^ is somewhat smaller than a typical rotation 
velocity for a galaxy fike the Milky Way {V,^^ ~ 220 km s"^), but 
this value is in line with expectations for the dark matter halo, once 
the effects of baryonic in-fall have been included [e.g., Klypin, Zhao, 
& Somerville 2002). 
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Fig. 12. — The average differential mass fraction, df /dMsz.t, nor- 
malized relative to host mass and satellite mass. The upper set 
of (bold) curves were computed for the n = 1 cosmology with 
A4ost = 1-4 X 10^2 at 2 = (solid) and Mhoat = lO" 

Mq (dotted), 3 X lO^^ Mq (long-dash) and 10^^ M© (dot-dash) 
all at z = 0.6. The lower set of thin curves correspond to BSI halos 
of Mhost = 1-4 X 10^2 M0 at 2 = (solid) and Mt^^st = 3 X lO^^ 
Mq at 2 = 0.6 (long-dash). The crosses reflect an analj^ic fit to the 
n = 1 results, as discussed in the text [see Eq. (14)]. 



massive than Mgat, /(< Afgat)- Figure 12 shows the differ- 
ential mass fraction df /dMsat, normalized relative to the 
host mass for several different host masses and redshifts 
in both the n = 1 and BSI models. The results are ap- 
proximately self-similar with respect to the host mass, and 
for the n = 1 case can be well-represented by the analytic 
form 

df 



dx 



X 


—a 


X 




exp 




_xo 




xo_ 



(14) 



with X = Msat/Mhost, a = 0.6 and xq = 0.07 ± 0.05. The 
quoted range in xq characterizes well the rms scatter from 
realization to realization (not shown) . This function (with 
Xq = 0.07) is shown as the set of crosses in Figure 12. As 
expected, the mass fractions are somewhat lower for the 
BSI model halos. The other CDM-type models all yield 
differential mass functions similar to those of the n = 1 
case. While in the next section we present results for a 
particular choice of host mass, the self-similarity demon- 
strated here implies that results at a fixed satellite-to-host 
mass ratio x, can be scaled in order to apply these results 
to any value of Mhost ■ 

4.3. Mass Fractions and Gravitational Lensing 

DKOl used flux ratios in multiply- imaged quasars to 
constrain the substructure content of galactic halos to be 
/ = 0.006 - 0.07 (at 90% confidence) for M^at^ 10*^ - 10^° 
Mq. In their sample of lens systems, the lens redshifts 
span the range 0.31^ Z£^ 0.97 with a median lens red- 
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Fig. 13. — The fraction of the parent halo mass that is bound up 
in substructure in the mass range between 10® Mq and Msat as a 
function of Msat • The host halo in each case has M = 3x10^2 Mq at 
z = 0.6. Lines reflect the mean over all realizations, and results 
are shown for the n = 1 model (solid), Rl model (dot-short-dash), 
iTg = 0.75 (dashed), as = 0.65 (dot-long-dash), and BSI (dotted). 
The error bars on the top set of lines reflect the 90 percentile range 
determined using 200 merger tree realizations for the n = 1 case 
(the other models in the top set of lines have very similar scatter). 
The bottom set of errors reflect the same range determined using 
50 realizations of the BSI model. 
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Fig. 14. — The fraction of mass in substructure in a central, 
cylindrical projection of radius 10 kpc, for the same set of halos 
as shown in Figure 13. The line types are the same as those in 
Figure 13, and again represent the mean fraction computed over all 
realizations. The large and small error bars represent the 90 and 64 
percentile ranges, respectively. A down-arrow is plotted instead of 
a lower, large error tick if at least 5% of the realizations had / = 
in that bin. A down-arrow with no accompanying lower error bar 
means that at least 18% of the realizations were without projected 
substructure in that bin. 



shift of Zi ~ 0.6. Our primary goal in this section is to 
make predictions aimed at lensing studies. Consequently, 
we present results for host systems at z = 0.6, and with 
A^host = 3 X 10^^ M0, which was taken as a typical lens 
mass in Dalai & Kochanek (2002, DK02). 

In Figure 13, we exhibit results for each of our inflation- 
derived power spectrum models. Here, /(lO^ Mq < M < 
Msat) is the cumulative fraction of host halo mass that 
is bound up in substructures with mass larger than 10^ 
Mq and less than Afgat- As expected from our discussion 
in §4.2, the mass fraction in substructure is not a strong 
function of the tilt of the primordial power spectrum, but 
it is sensitive to a sudden break in power at small scales. 
Specifically, the subhalo mass fraction in the BSI model is 
roughly a factor of ~ 3 below that seen for the CDM-type 
spectra in this mass range. The top set of error bars reflect 
the 90 percentile range derived using 200 realizations for 
the n = 1 model (other CDM-type models show similar 
scatter) and the bottom set of errors reflect the same range 
determined from 50 realizations of the BSI spectrum. 

Rather than the total mass fraction, lensing measure- 
ments are sensitive to the mass fraction in substructure 
projected onto the plane of the lens at a halo-centric dis- 
tance of order the Einstein radius of the lens. Re '--^ 5 — 15 
kpc. In Figure 14 we plot /sat(> 10^ M©) projected 
through a cylinder of radius 10 kpc centered on the host 
halo for the same set of halos shown in Figure 13. The 



large and small error bars reflect the 90 and 64 percentile 
ranges, respectively, in measured projected mass fractions 
derived using 200 n = 1 realizations (top set) and 50 BSI 
realizations (bottom set). A down-arrow is plotted instead 
of a lower, large error tick if at least 5% of the realizations 
had / = in that bin. A down-arrow with no accompa- 
nying lower error bar indicates that at least 18% of the 
realizations were without projected substructure in that 
bin. 

The projected mass fractions are not as severely sup- 
pressed relative to the volume-averaged mass fractions as 
one might expect given that tidal forces act systematically 
to destroy substructure near host halo centers (see Figure 
5). The reason is that we are examining substructure in 
a cylindrical volume, and picking up subhalos with large 
halo-centric radii. We illustrate this effect in Figure 15, 
where we compare the mass fraction in cylindrical projec- 
tion radius p with the mass fraction in spherical shells with 
the same value of spherical radius r. Notice that the mass 
fraction in spherical regions is significantly reduced in the 
center, while the projected mass fraction is less severely af- 
fected. Of course, the mass fraction approaches the global 
value at large radii. Figures 16 and 17 demonstrate how 
the mass fractions change as a function of projection ra- 
dius for various subhalo mass cuts for the n = 1 and BSI 
models respectively. Notice that the relative drop in mass 
fraction as a function of projection radius is more pro- 
nounced in the BSI model than in the n = 1 case. This 
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Fig. 15. — The cumulative mass fraction of substructure with 10® 
M0 < Msat < 10^ M0 shown spherically averaged as a function of 
radius r (solid line) and in projection, as a function of the projected 
radius p (dashed line). The averages (lines) and 64% percent ranges 
(error-bars) were determined using 200 realizations of M = 3 X lO^'^ 
Mq halos z = 0.6 for our n = 1 model. Down-arrows indicate that 
more than 18% of the realizations had / = in the corresponding 
radial bin. 
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Fig. 16. — Mass fraction in substructure in cylindrical projection 
of radius p for the same set of n = 1 halos described in Figure 15. 
The upper left, upper right, and lower left panels show the mass 
fraction profiles in subhalos larger than 10® Mq and less than 10^, 
10* and lO'^ Mq respectively. Error bars reflect the same percentile 
ranges as do those in Figure 15. The bottom right panel illustrates 
how the mean projected mass fraction profiles vary as a function of 
the maximum subhalo mass considered: 10®-^, 10^, 10*, lO'^, and 
10^° Mq from bottom to top. 



reflects the fact that tidal disruption is more important in 
the BSI case and core-like behavior of the subhalo radial 
distribution sets in at a larger radius in this model. 

4.4. Warm Dark Matter and Gravitational Lensing 

In the previous section we demonstrated that the sub- 
structure mass fraction is sensitive to abrupt changes in 
the power spectrum and used the BSI model as a specific 
example. In this section we investigate these differences in 
the context of WDM. We label the different WDM models 
by the warm particle mass mw, and assume the canonical 
case of a "neutrino-like" thermal relic with two internal 
degrees of freedom, = 2. 

Figure 18 shows the total mass fraction of 3 x 10^^ M0 
host halos at z = 0.6 as a function of Afgat implied by 
our three WDM model power spectra compared to our 
standard ACDM case. For substructure smaller than ^ 
10^ Mq, the differences between the models are as large as 
an order of magnitude or more, and even the largest WDM 
particle mass (3 keV) provides a potentially measurable 
suppression of substructure. Figure 19 shows the mass 
fraction in projected cylinders of radius 10 kpc. 

The differences in mass fractions seen for the different 
models in Figures 18 and 19 come about because subhalos 
become less concentrated relative to their host halos as the 
WDM particle mass is decreased and power is suppressed 
on larger scales, much like the BSI case. In true WDM 
models there are additional processes that, in principle, 
can alter the formation and density structure of dark mat- 
ter halos. In Figures 18 and 19, we have only accounted 



for the effect of the power spectrum on substructure mass 
fractions and assumed that the density structure of WDM 
halos is identical to that for CDM halos. For high mass sys- 
tems, this is a sensible approximation (Colin et al. 2000a; 
Avila-Reese et al. 2001), but this approximation should 
break down at small masses and lead to further suppres- 
sion of substructure. 

One consequence of a WDM particle with non-negligible 
velocity dispersion is that gravitational clustering is re- 
sisted by structures below the effective Jeans mass of the 
warm particles {e.g., Hogan & Dalcanton 2000; Bode et al. 
2001): 

(15) 

For both the m-yv = 1-5 keV and mw — 3.0 keV models, 
Mj <C IO^Mq when 10, so all halos of interest in this 
context are minimally affected. The situation is some- 
what more complicated in the m^v = 0.75 keV model, 
where AfjJ^ 10^ Mq for redshifts 2. We therefore ex- 
pect that the formation of these halos should be suppressed 
compared to the predictions of the EPS formalism. This 
suppression should only have a minor effect on our pre- 
dictions because we restrict ourselves to satellite masses 
10^ M0 and most surviving subhalos are accreted at 
2. In the interest of simplicity, we chose to ignore this 
effect here. As a result, we may significantly over-predict 
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Fig. 17. — Mass fraction in substructure in cylindrical projection 
for the BSI model. All line types and panels are as described in 
Figure 16. 



substructure mass fractions at low Mgat in these cases. In 
the context of this study, this is a conservative approach 
because the true mass fraction would be reduced by these 
effects, bringing it further away from the measured sub- 
structure mass fractions and standard ACDM predictions. 

In addition to the effective Jeans suppression, WDM 
halos, unlike their CDM counterparts, cannot achieve ex- 
tremely high densities in their centers due to phase space 
constraints (Tremaine & Gunn 1979). In the early Uni- 
verse the primordial phase space distribution of the WDM 
particles is a Fermi-Dirac distribution with a maximum of 
/max = ffw/^pL low energies (/ipL is Planck's constant). 
For a coUisionless species, the phase space density is con- 
served and this maximum phase space density may not be 
exceeded within WDM halos. If we define the phase den- 
sity as Q = p/{2na^)^/^ then the maximum allowed phase 
density is (Hogan & Dalcanton 2000) 



5.2 X 10" 




Mo/pc3 



(16) 



(km/s)3 

This limit implies that WDM halos cannot achieve the cen- 
tral density cusps of the kind observed in simulated CDM 
halos. Instead, we expect a core in the density profile. For 
viable WDM models, the phase space core is expected to 
be dynamically unimportant for any halo massive enough 
to host a visible galaxy (ABW). However, for the lowest- 
mass subhalos (M^ 10^ ^q) the presence of phase space- 
limited cores may be important because halos with large 
cores are less resistant to tidal forces than cuspy halos. 

We have attempted to estimate (crudely) how the phase 
space limit affects the substructure population of WDM 
halos by adopting our standard model of halo accretion 
and orbital evolution, but allowing the density structure 
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Fig. 18. — Total cumulative mass fractions in substructure more 
massive than 10^ Mq for the n = 1 and 3 WDM models. The 
models are ACDM (solid), mw = 3.0 keV (dashed), mw = 1-5 keV 
(dash-dot), and = 0.75 keV (dotted). For clarity, we show error 
bars only for the ACDM and mw = 1-5 keV models. The error bars 
and arrows have the same meaning as in Figure 14. 



of the appropriately small subhalos to be set by the phase 
space limit. For these calculations we used the phenomeno- 
logical density profile of Burkert (1995), 



PB{r) 



Po 



[l + r/rB][l + {r/rBy 



(17) 



The Burkert profile resembles the NFW form at large ra- 
dius, but features a constant density core at its center, 
and thus a velocity dispersion that approaches a constant 
at small r: ~ 0.55T4iax-^ For Burkert profiles, l^^x — 
0.86V^^ij.CB/gB(cB), where cb = Rvii/fB is the Burkert con- 
centration and (7b (y) = ln(l-|-?/^) -|-2 ln(l-|- y) — 2 tan~^(y). 
Solving for the phase density in the core (r ^ tb) and 
equating it with the maximum phase density of equation 
(16) yields the following relation for the maximum attain- 
able value of Cb: 




We assigned Burkert concentrations according to the fol- 
lowing prescription. First, we computed NFW concentra- 
tions Cvir, for each halo according to the BOl model. We 
converted from NFW concentration to Burkert concentra- 
tion Cb, by interpreting the BOl value of as the ra- 
dius at which dln/3(r)/dlnr|.r=rs = ~2. This implies that 

^ The numerical coefficient in Eq. (17) of ABW should be « 0.3 
rather than 0.2. 
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Fig. 19. — Cumulative mass fractions in substructure more mas- 
sive than 10^ Mq within a projected radius of 10 kpc for the same 
models shown in Figure 18. The linetypes are the same as in Figure 
18. 
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Fig. 20. — Total cumulative mass fractions in substructure more 
massive than 10^ Mq both with (dashed) and without (solid) esti- 
mating the effects of the phase space limit. The upper set of lines 
corresponds to the myv = 1.5 keV model and the lower set of lines 
corresponds to the mw = 0.75 keV model. The error bars and 
arrows are as in Figure 14. 



tb — 0.66rs or cb — l.Scvir- With this correspondence, the 
adopted Burkert profile achieves the maximum of its rota- 
tion curve at rj^ax — 0-99''max, where rmax is the radius at 
which the corresponding NFW halo achieves Knax- Simi- 
larly, Knax of the adopted Burkert profile is within 10% of 
the corresponding NFW Vmax for all relevant concentra- 
tions (1 < Cvir < 25). Second, we computed the maximum 
value of Cb allowed by the phase space constraints using 
Eq. (18). We then assigned each halo the smaller of these 
two values of Cb at the time of accretion. In this way, we 
guaranteed that the phase space constraint was met by all 
halos. We have checked that this prescription for Burk- 
ert halos does not yield any systematic bias in our results 
by applying it all of our CDM models. We found that it 
gave nearly identical results to that of our standard NFW 
model, which is not surprising in the context of our model 
and disruption criteria. 

We present our estimates of cumulative mass fractions 
in WDM models, including the effect of the phase space 
constraint, in Figure 20. It is clear, at least from this rough 
estimate, that the Tremaine-Gunn limit plays an impor- 
tant role only for the most extreme WDM models, 1 
keV, and only the smallest halos, ^ 10® M0. However, we 
emphasize that our new assumptions about WDM halos 
have not been tested with N-body simulations. Simula- 
tions have yet to examine the density structure of halos 
that saturate the phase space bound and most studies 
have ignored the initial velocity dispersion of the WDM 
particles (Colfn et al. 2000a; Avila-Reese et al. 2001, 
Knebe et al. 2002), but the Burkert profile assumption 
seems plausible. With these precautions in mind. Figure 
18 may be regarded as an approximate upper-limit on the 
substructure mass fraction for WDM halos. Any phase 



space bound or the effects of primordial velocity disper- 
sions on halo formation and density structure should lead 
to enhanced disruption, resulting in lower mass fractions. 

One physical process that might affect WDM (and BSI) 
models that we have not considered is top-down fragmen- 
tation {e.g., Knebe et al. 2002). It is possible that power 
can be transported from large scales to small in trun- 
cated models, resulting in a population of low-mass halos 
that would not be accounted for in Press-Schechter theory. 
While such a process could result in a higher substructure 
abundance than that estimated using our model, there are 
reasons to believe that the effect should be fairly small. 
Systems that form in this manner collapse quite late, and 
their density structure likely would be very diffuse com- 
pared to their hierarchically-formed brethren. Therefore, 
it is less likely that systems formed via fragmentation could 
survive tidal disruption once incorporated into a galactic 
halo. 

4.5. The Dwarf Satellite Problem 

Comparisons between the predicted subhalo population 
and the observed dwarf galaxy abundance are usually made 
by comparing counts as a function of maximum circular 
velocity, Vmax- This is a sensible mode of comparison be- 
cause it sidesteps the complicated issues of star formation 
and feedback in these poorly-understood galaxies. Yet, 
there are considerable uncertainties, even for this method 
of comparison, and it is likely that efforts to compare pre- 
dictions as a function of dwarf luminosity (Somerville 2002; 
Benson et al. 2002) in tandem with velocity comparisons 
will be needed in order to fully understand the nature of 
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Fig. 21. — The lower and upper sets of points show a scatter 
plot of Vmax and rmax for all surviving halos produced in 10 merger 
tree realizations for the n = I, and low-normalization, ag = 0.65 
models respectively. The thick solid line shows the locus of points 
in the Vmax-r'max plane that corresponds to the central value of 
the measured velocity dispersion of Carina (cr^ = 6.8 ± 1.6 km s~^) 
given the measured King profile parameters of Carina (Table 2) 
and assuming NFW halos. The thin solid lines correspond to the 
ztlcr quoted errors in velocity dispersion for Carina. The dashed 
line corresponds to the locus of points that is consistent with the 
central value of the measured velocity dispersion of Draco (cr* = 
9.5 km s^i). 



the dwarf sateUite problem. 

For most satellites, the quantity that is observed and 
used to infer the halo Knax is the line-of-sight stellar ve- 
locity dispersion, ct*. As discussed in S02 and H03, the 
mapping between cr* and T^nax depends upon the theo- 
retical expectation for the density profile of the subhalo 
as well as on the stellar mass distribution of the galaxy. 
An additional complication concerns the unknown velocity 
anisotropy of the stars in the system. 

A phenomenologically-motivated approximation for the 
stellar distribution in a dwarf galaxy is the spherically 
symmetric King profile (King 1962), 



(19) 



where z = [1 + (r/rc)^]/[l + (rt/rc)^], and rt are the 
core and tidal radii of the King profile, and /0*(r > rt) = 0. 
The normalization is not important in what follows. 

If we assume that a stellar system described by Equa- 
tion (19) is in equilibrium and embedded in a spherically 
symmetric dark matter potential characterized by the cir- 
cular velocity profile Vc{r), then the radial stellar velocity 
dispersion profile crrir), can be computed via the Jeans 
equation: 

d[/0*cr2] 



where the anisotropy parameter (3=1 — a\/2a'^, a± is the 
tangential velocity dispersion, and /3 = corresponds to 
an isotropic dispersion tensor. A measured, line-of-sight 
velocity dispersion is determined by the projected velocity 
dispersion profile weighted by the luminosity distribution 
sampled along the line-of-sight. For the isotropic case it is 
given by 

)K^(r')dr' ^ ^ 

) c\ J ^21) 



/o'V*(?'')dr' 

assuming a constant mass-to-light ratio. If a galaxy has 
a measured stellar profile (the King parameters in this 
case) and measured value of cr^,, then the Jeans equation 
places only one constraint on the rotation curve of the 
system, Vc{r). We expect the halo velocity profile to be 
at least a two-parameter function {e.g., the NFW profile) 
so determining Vmax requires some theoretical input for 
the expected form of Vc{r) in order to provide a second 
constraint. 

Motivated by dark matter models, we assume that the 
global rotation curve is set by an NFW profile associated 
with the dwarf galaxy halo. The rotation curve for an 
NFW halo is fully described by specifying two parameters 
and a natural pair is T^ax and rmax- For any given cos- 
mology, the relation between Vmax and rmax is expected to 
be rather tight, and this provides a second, theoretically- 
motivated constraint that sets the Knax-f* mapping im- 
phed by Eq. (20). 

The Vmax-rmax relationships for surviving subhalos in 
two of our models are shown in Figure 21. The lower set 
of points corresponds to our standard n — 1 model and 
the higher set of points is derived from our ug = 0.65 
model. In each case, we plot one point for each surviv- 
ing halo in 10 model realizations. The strong correlation, 
^max c>c V2axi7 — 1-3, foUows directly from the input cor- 
relations between Mvir(z), and Cvir (see §2.3 and BOl). 
The normalizations and slopes are influenced by the cos- 
mology, accretion times and (mildly) by the orbital history 
of the subhalos.^" 

The thick solid and dashed lines in Figure 21 show the 
locus of points in the Vmax-rmax plane that correspond 
to the central values of the observed line-of-sight velocity 
dispersions for Carina and Draco respectively, given their 
measured King profile parameters. Our adopted cr* val- 
ues and King parameters are listed in Table 2 along with 
appropriate references. The light solid lines illustrate how 
these contours expand when we include the ztlcr measure- 
ment error in cr^ for Carina. A similar (although narrower) 
band exists for Draco, but we have omitted it for the sake 
of clarity. Consistency with the observed King parameters 
and velocity dispersions requires each dwarf to reside in 
a halo with structural parameters that lie within the re- 
gion of overlap between the contours and the model points. 
For example, in the n — 1 model Carina is expected to 



reside in a halo with V,^ 



11 km s and r„ 



1 



kpc. For the as = 0.65 model, Carina is expected to sit 



in a larger halo, with Vm 



dr 



-p.{r)V,'ir) - 2P{r)p,ir)a^{r), (20) 



29 km s ^ and r^ax ~ 10 
kpc. Similar comparisons hold for Draco and all of the 

The scatter in the V'max-f'max plane should be larger than that 
shown here because we have not included the expected scatter in 
the input Cvir-Mvir relation. For CT(logCvir) ~ 0.14 (BOl; W02), the 
implied scatter is a'(logrinax) ~ 0.18 at fixed Vmax- For (T(logCvir) ~ 
0.08 (Jing 2000), the implication is (T{logrmax) ~ 0.11. 
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Local Group dwarf satellites and these comparisons can 
bo made in a similar way for any cosmology. The point is 
that the maximum velocities that are assigned to satellite 
galaxies are cosmology-dependent. Therefore, "observed" 
velocity functions are also cosmology-dependent because 
theoretical inputs are used to convert from cr* to Vmax- 

In Table 2 we show estimates for halo Knax values for 
the observed Milky Way satellites under the assumption 
that /? = along with a similar analysis for /? = 0.15 (val- 
ues in parentheses). We estimate halo Knax values for six 
different power spectra, relying on the model-dependent 
''max-Knax relationship for substructure in each case, and 
taking the central values of the measured velocity disper- 
sions for each halo. Taking the quoted ±la range for the 
measured velocity dispersions typically leads to a shift in 
Vmax of ~ 30% which is considerable compared to the in- 
herent scatter in the Myjr — Cvir relation. For reference, we 
have also included the adopted Knax values from the origi- 
nal K99 work on the dwarf satellite problem. As expected, 
the implied Knax values become larger as we explore mod- 
els with less galactic-scale power. Our estimates for the 
n = 1, /3 = case are close to those of K99. 

The left panel of Figure 22 shows the Milky Way satel- 
lite counts for each model, assuming = 0, along with 
the predicted velocity functions for each model. In addi- 
tion to the satellites listed in upper portion of Table 2, 
we have also included the Small Magellanic Cloud (SMC, 
with Knax = 60 km s~^, estimated by Stanimirovic 2000 
to include a substantial contribution from the baryonic 
component) and the Large Magellanic Cloud (LMC, with 
Knax = 50 km s^^, van der Marel et al. 2002) in our cu- 
mulative velocity functions. In the standard case (n = 1) 
the discrepancy sets in at I^ax ~ 30 km s~^, requiring 
roughly one-in-ten halos with Knax ^ 10 — 20 km s~ 
to be the host of an observed galaxy. The extremely 
tilted model with cjg ~ 0.65 actually under-predicts the 
dwarf count for large systems. Interestingly, dwarfs in the 
dn/dlnfc = —0.03 RI model as well as the ag = 0.75 case 
are consistent with inhabiting the ten most massive sub- 
halos, with only Sextans standing as an outlier. The BSI 
model also looks to be in good agreement with the data 
for VJnax^ 12 km s""'^, lending support to the conclusions 
of Kamionkowski & Liddle (2000). However, the prob- 
lem of Local Group satellites is not completely "solved" 
in any of these models because the velocity function con- 
tinues to rise below the velocity-scale of Sextans in all 
cases. What changes in the low-power models is the na- 
ture of the discrepancy. In one extreme, the mismatch 
sets in at Knax ~ 30 km s~^ and gradually becomes worse 
for smaller systems. In the other extreme, the mismatch 
seems to imply a sharp threshold for dwarf galaxy forma- 
tion at a scale near T^nax ^10 — 20 km s^^. 

Unfortunately, a detailed accounting of the mismatch 
is difficult, even for a given cosmology. The dwarf Knax 
estimate is very sensitive to the velocity anisotropy param- 
eter, /3. For example, assuming /? = 0.15 leads to VJnax val- 
ues that are significantly lower than in the isotropic case 
because rotational support has been traded for pressure 
support. The velocity function comparisons with /? = 0.15 
are shown in the right panel of Figure 22 and the Knax val- 
ues are listed in parentheses in Table 2. In this case, only 
the (Tg — 0.65 case and the BSI model can account for the 



dwarf population without a differential feedback mecha- 
nism. The rest of the models over-predict the counts, with 
the greatest apparent discord in the n = 1 case. The spe- 
cific choice of /? = 0.15 serves mainly to illustrate the effect 
of a minor anisotropy. We chose this value because it is 
typical of what is seen in the central regions of simulated 
dark matter halos (Colin et al. 2000b), and therefore it 
seems a reasonable possibility for the anisotropy parame- 
ter of particles in dwarf galaxies. 

5. CAVEATS 

In this study we employed a simple, semi-analytic model 
based on many previous studies (LC93; SK99; BKW; TBOl; 
BOl; W02; IIFM02) and designed to produce large num- 
bers of halo realizations with minimal computational ef- 
fort. In developing this model, we have made many sim- 
plifying assumptions. In this section we draw attention to 
many of these shortcomings and discuss how they might 
affect our results and be improved upon in future work. 

Among the most obvious omissions in this work is the 
neglect of any disk or bulge component in each halo. We 
have specifically chosen to ignore the effects of central 
galaxies because the physics of dark halo formation is 
relatively well-understood compared with that of galaxy 
formation. This allows us to ground our work against 
dissipationless N-body simulations. In order to include 
a galactic component, one is forced to adopt many poorly- 
constrained models and assumptions regarding gas accre- 
tion, cooling, angular momentum distributions, feedback, 
and the effects of substructure on the host galaxy itself. 
Once a reliable framework for the dark matter has been 
developed, we can use this as a foundation for more specu- 
lative (yet interesting) explorations involving the baryonic 
components. 

A central (disk) galaxy would add to the dynamical fric- 
tion force experienced by sublialos orbiting near the plane 
of the disk and cause halos on highly inclined orbits to be 
tidally heated during rapid encounters with the disk po- 
tential (e.g., Gnedin & Ostriker 1999; Gnedin, Hernquist, 
& Ostriker 1999; TBOl). These effects lead to enhanced 
satellite disruption. Conversely, subsystems that are mas- 
sive enough to host galaxies might be rendered more re- 
sistant to tidal disruption because their central densities 
would be enhanced by the presence of cool baryons. For 
low-mass halos, the net effect of a central galaxy would 
likely be to reduce the substructure count, mainly at small 
radii. Even without including these effects, we find that 
the substructure fraction drops significantly at small radii 
because of the dark matter potential, and that a large part 
of the projected central mass fraction comes from subhalos 
at large radii that are picked up in projection. Neverthe- 
less, projected mass fractions are rather sensitive to the 
size of the core in the subhalo radial distribution (Chen et 
al. 2003), so if the core region were larger as a result of 
a central galaxy, the implied lensing signal would be re- 
duced relative to our estimates. We find that eliminating 
all substructure within 20 kpc of the halo center, reduces 
the projected mass fractions in subhalos less massive than 
lOio Mq, /lo, by ~ 30%. 
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Fig. 22. — Satellite halo velocity functions for six of our models compared with the velocity functions of the Milky Way satellites after 
accounting for the cosmology-dependent mapping between cr* and Vmax. The squares represent the velocity function of Milky Way satellites 
based on the data in Table 2. The lines represent the means over all realizations and the error bars reflect the dispersion among these 
realizations. The models are labeled in each panel. 



Table 2 

Characteristics of Milky Way satellites. 



Satellite 




Tc 


rt 


Vmax 


Vmax 


Vmax 


Vinax 


Vmax 


Vmax 


Vmax 










K99 


n=l 


o-g = 0.83 


RI model 


(Tg = -75 


BSI 


(T8 = 0.65 




( km s-1) 


(kpc) 


(kpc) 


( km s-l) 


( km s-i) 


( km s-i) 


( km s-i) 


( km s-1) 


( km s-l) 


( km s-l) 


Sagittarius 


11.7 ±0.7 


0.44 


3.0 


20 


17 (13) 


20 (15) 


23 (17) 


24 (17) 


26 (19) 


35 (25) 


Fornax 


10.5 ± 1.5 


0.46 


2.3 


18 


15 (11) 


18 (13) 


20 (15) 


21 (15) 


23 (17) 


30 (21) 


Draco 


9.5 ± 1.6 


0.18 


0.93 


17 


21 (15) 


31 (20) 


38 (25) 


41 (26) 


41 (29) 


67 (44) 


Ursa Minor 


9.3 ± 1.8 


0.20 


0.64 


16 


21 (15) 


31 (20) 


37 (25) 


41 (26) 


42 (29) 


69 (43) 


Leo I 


8.8 ±0.9 


0.22 


0.82 


15 


17 (12) 


24 (16) 


28 (19) 


30 (20) 


33 (23) 


51 (32) 


Carina 


6.8 ± 1.6 


0.21 


0.69 


12 


11 (8) 


14 (10) 


18 (13) 


18 (13) 


22 (16) 


29 (19) 


Leo II 


6.7 ± 1.1 


0.16 


0.48 


12 


13 (9) 


19 (13) 


23 (16) 


24 (16) 


27 (19) 


40 (26) 


Sculptor 


6.6 ±0.7 


0.11 


1.5 


11 


13 (9) 


18 (12) 


21 (15) 


23 (15) 


26 (18) 


38 (24) 


Sextans 


6.6 ±0.7 


0.34 


3.2 


11 


8 (7) 


10 (7) 


10 (8) 


11 (8) 


12 (9) 


14 (10) 


SMC 








30 


< 60* 


< 60* 


< 60* 


< 60* 


< 60* 


< 60* 


LMC 








50 


50 


50 


50 


50 


50* 


50* 



Note. — The names of the satellite galaxies of the Milky Way are given in column (1). We consider only those galaxies with galactocentric 
distances smaller than 300 kpc. Columns (2-4) give the measured line-of-sight velocity dispersion, and King core and tidal radii for each satellite. 
The exceptions are the LMC and the SMC. These galaxies have measured rotation speeds, listed in Columns (6-11). All velocities are expressed 
in units of km and distances are in kpc. Column (5) gives the value of Vmax assigned to the halo of each satellite by by K99. Columns 
(6-11) and above the horizontal space give the value of Vmax that we estimated for each satellite based on its measured velocity dispersion and 
King profile parameters and assuming the primordial power spectra specified at the top of each column. Values of Vmax listed without parentheses 
were calculated assuming an anisotropy parameter of /3 = 0, and those inside parentheses assume /3 = 0.15. Except for the case of Draco, we use 
the velocity dispersions and King profile core and tidal radii for the Milky Way satellites given in the review article by Mateo (1998). For Draco, 
we use the parameters quoted by Odenkirchen et al. (2001). The quoted maximum rotation speeds for the LMC and SMC were taken from van 
der Marel et al. (2002) and Stanimirovic (2000). The LMC rotation curve is observed to be flat from 4kpc out to > 8.9kpc. For many of the 
low-power models (indicated by superscript "*"), the flat portion of the curve (~ rmax) is expected to be at larger radius. In order to explain this 
in the context of these models, we must suppose that baryonic in- fall plays an important role in setting the properties of dark matter rotation 
curves (Blumenthal et al. 1986). In this case, the measured value of Vmax (j'max) is larger (smaller) than it would be for the pristine halo prior to 
baryonic contraction. The SMC rotation curve is even more likely to be influenced by baryons (Stanimirovic 2000), and baryonic in-fall is likely to 
be of some importance for all cases. While not demanded by the data, the effects of baryonic in-fall could be important for all satellites, thus the 
listed Vmax values should be considered lower limits. Lastly, the large value of rmax associated with Draco in the trg = 0.65 case may be difficult 
to reconcile with the kinematic data of Kleyna et al. (2002) and may disfavor a model with such low power. 
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In addition to the considerable uncertainties associated 
with galaxy formation, there are potential shortcomings in 
our efforts to model substructure properties in the context 
of collisionless dark matter physics alone. For example, we 
have allowed for only a mild redistribution of mass within 
the tidal radii of the orbiting subhalos up until the time the 
subhalo is totally disrupted. The work of S02 and H03 sug- 
gests that this effect may be larger, but these results may 
have been compromised by limited numerical resolution or 
inappropriate assumptions regarding initial subhalo orbits 
and/or accretion times. In this sense, our approach repre- 
sents a conservative extreme because we assume that the 
surviving subhalo density structure is typically very simi- 
lar to that of halos in the field. We study the effects of tidal 
mass redistribution in a forthcoming paper (.1. Bullock, K. 
Johnston, & A. Zentner, in preparation). In addition, we 
have adopted a halo concentration relation (BOl) that has 
not been confirmed for 10*^ Mq. Similarly, the EPS 
merger tree calculations have yet to be tested in the very 
low-mass regime. In light of these extrapolations, it is im- 
perative that our results be tested, and updated using the 
next generation of numerical studies. 

Finally, our model does not treat the substructure pop- 
ulation self-consistently. We have neglected any subhalo- 
subhalo interactions which could serve to increase the in- 
ternal heating of substructure and modify dynamical fric- 
tion timescales as orbital energy is exchanged between sub- 
halos and traded for internal energy. We have adopted 
the approximation that all in-falling halos are "distinct" 
and have no subhalos of their own (see Taylor & Babul 
2003 for a study of merger tree "pruning"). However, 
our "tree-level" calculations suggest that the substructure 
mass fraction is uniformly ~ 10% regardless of host mass, 
so we expect the that this correction would typically af- 
fect our derived mass fractions by ^ 10%. Considering the 
assumptions that have gone into our calculations and the 
current level of observational precision, this level of error 
is acceptable; however, it may need to be improved upon 
as observations zero in on the masses of the subclumps 
responsible for the lensing signals and the mass fractions 
in these subclumps. 

6. CONCLUSIONS AND DISCUSSION 

The abundance of substructure in dark matter halos is 

determined by a continuous competition between accre- 
tion and disruption. Accreted subhalos with dense cores 
are resistant to disruption, but over time their orbits de- 
cay, their mass is stripped away, and they are often de- 
stroyed. The model we presented here allows us to follow 
the complicated interplay between density structure, or- 
bital evolution, and accretion time in order to determine 
how changes in the power spectrum affect the final sub- 
structure population in galaxy-sized halos. For a fixed set 
of cosmological parameters, changes in the power spectrum 
manifest themselves by changing collapse times for halos, 
where less power leads to later accretion times and lower 
densities. We have specifically focused on tilted models 
that help to relieve the central density crisis facing CDM 
and that may be favored by joint CMB and large-scale 
structure analyses (Spergel et al. 2003). We have also 
considered a BSI inflation model and WDM models, where 
the power is sharply reduced on small scales. 



For a large class of CDM-type models, including models 
with significant tilt and running, we find that the fraction 
of mass bound up in substructure /, for galaxy- mass ha- 
los is relatively insensitive to the slope of the primordial 
power spectrum. This is because both the host halos and 
their accreted subsystems collapse later in these models, 
in a roughly self-similar way as power is reduced. Note 
that this result would have been roughly expected if we 
were varying only the overall normalization in the mod- 
els because the relative redshifts of collapse would be in- 
variant (we assume host halos are small enough that they 
collapse before ^; ^ 0). Our investigation suggests that 
this intuitive description holds even for tilted and running 
index models, at least over the parameter range we have 
explored. All indications are that this insensitivity to the 
tilt of the power spectrum is a rather robust result, and 
should hold even if some unknown factor has caused our 
overall normalization in predicted mass fractions to be in 
error {e.g., our exclusion of central galaxies). Interest- 
ingly, the shape of the mass function, f{x = Afgat/Mhost), 
is also relatively insensitive to the mass of the host halo 
[Eq. (14)], and a similar shape holds for all of the tilted 
models we explored. 

The similarity in mass fractions breaks down for models 
with sharp features in their power spectra, like our BSI and 
WDM models. In these models, low-mass halo formation is 
delayed significantly relative to the formation time of their 
hosts. Consequently, fewer subhalos are dense enough to 
withstand the tidal field they experience upon accretion. 
We find that for the relevant WDM and BSI models, the 
mass fraction in substructure is reduced by a factor of ^ 3 
compared to the standard/tilted ACDM models. 

Inspired by recent attempts to measure substructure 
mass fractions using multiply-imaged quasars, we applied 
our model to ensembles of host halos with M = 3 x 10^^ 
M0 at z = 0.6, which represents the expectation for mas- 
sive lens galaxies (DKOl; DK02). For the ACDM/tilted 
cases, we found substructure mass fractions within a 10 
kpc projected radius in systems less massive than M = 
10*, 10^ and IQ^^Mq of /g ~ 0.2 - 0.4%, /g ~ 0.4 - 
1.5%, and /lo — 0.6 — 2.5% at the 64 percentile range. 
These estimates are consistent with, but on the low side 
of, first attempts to measure the substructure fraction us- 
ing multiply- imaged quasars by DKOl, who obtain / ~ 
0.6% — 7% at 90% confidence, with an upper-mass limit 
of 10* - 10^° Mq (N. Dalai, private communication). The 
lensing results disfavor the BSI model, which leads to mass 
fractions /s =± 0.01 - 0.06%, /g ~ 0.02 - 0.2%, and fw ~ 
0.03 — 0.4% at 64%. This is true unless the break scale in 
the power spectrum is pushed to such a small value that 
this model no longer has the attractive feature of alleviat- 
ing the central density problem. A tow = 0.75 keV WDM 
model is similarly disfavored, and even our highest mass 
WDM case, mw = 3 keV, has a typical projected fraction 
(/g ~ 0.4%) that is low compared to the DKOl estimate. 
Again, this indicates that if the warm particle is a thermal 
relic, the mass must be large enough that WDM no longer 
mitigates the small-scale problems of standard CDM. Yet, 
these results are interesting because they show how lensing 

Note that these self-similar trends brealc down on cluster-mass 
scales since recent accretion, which is a strong function of the overall 
power, likely plays an important role in these objects. 
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may be used as one of the few probes of the WDM particle 
mass in the range ^ 1 keV or a break in the primordial 
power spectrum at large wavenumber. 

Clearly these conclusions must be regarded with some 
caution. In addition to the uncertainties of modeling dis- 
cussed in §5, other issues make drawing definite conclu- 
sions difficult. For example, we have only accounted for 
the substructure within the virial radius of the host halo, 
yet the anomalous flux ratios of lensed images are sensi- 
tive to the presence of small halos along the line-of-sight 
to the lens. Keeton (2003) showed that field halos can 
have a significant lensing effect even if they are separated 
from the lens by several tenths in rcdshift and in hierar- 
chical, CDM-type models, small field halos arc ubiquitous. 
Chen et al. (2003) showed that the relative effect from ha- 
los outside the virial radius of the lens is typically a few 
percent, but may be as large as 20 — 30% of that from 
subhalos, depending upon assumptions about the subhalo 
population. Also, as the mass fraction in substructure of 
a given mass depends on the mass of the host [Eq. (14)], 
it may be important to constrain the host halo mass in or- 
der to fully exploit the ability of lensing measurements of 
substructure to probe cosmology and structure formation. 

We compared our model predictions for the cumula- 
tive subhalo velocity function, A''(> X^nax), to the satellite 
galaxy count of the Milky Way. The approach here was to 
estimate Knax for each satellite galaxy's dark matter halo 
based on its observed line-of-sight velocity dispersion, tr*. 
We emphasized in §4.5 that the mapping between cr* and 
T4iax is sensitive to theoretical prejudice regarding the den- 
sity structure of the dwarf galaxy's halo as well as the un- 
known velocity anisotropy parameter of the system, p. For 
a fixed value of /?, less concentrated host halos imply larger 
values of Knax because halo rotation curves are more slowly 
rising and stars probe only the inner ~ 1 kpc of the halo. 
Interestingly, this implies that tilted models and truncated 
models, do significantly better than n = 1, ACDM in re- 
producing apparent dwarf counts, even though their mass 
fractions are similar. While our estimates of Vmax cannot 
be considered robust because of the simplicity of our model 
and the fact that the cr^ — V,nax mapping is very sensitive 
to the inner structure of the subhalos, the general trends 
that we illustrate should persist in more elaborate stud- 
ies and N-body simulations. Moreover, our resTilts reveal 
yet another reason why it is difficult to consider the dwarf 
satellite problem a serious challenge to CDM theory, the 
nature of the dwarf satellite problem is very sensitive to 
cosmology, the power spectrum, and assumptions about 
the shape of the velocity ellipsoid for stars in dwarf galax- 
ies. 

When we fix /? = 0, our dn/dlnk = —0.03 RI model, 
0-8 = 0.75 model, and the BSI case all do well in match- 
ing the known satellite population of the Milky Way for 
Knax^ 20 km s^ . Our lowest power model (o"8 = 0.65, 
n ~ 0.84, and mild running) actually under-predicts the 
dwarf count for Knax~ 30 km s~^. However, this result 
is achieved only for the optimistic assumption of isotropic 
velocities. If we adopt a small level of anisotropy, (i = 0.15, 
consistent with the centers of simulated dark matter ha- 
los, agreement for most models is worsened. Only the BSI 
and (Tg = 0.65 models show good agreement in this case. 
Yet, even with /3 = 0.15, the RI and o-g = 0.75 models still 



compare more favorably than the n = 1 case with /3 = 0. 

What do these results imply for the dwarf satellite prob- 
lem? In all models, including those with truncated power, 
the velocity function of subhalos continues to rise below 
the scale of the smallest observed Milky Way satellite, 
^max~ 10 km s~^. No matter how one modifies the power 
spectrum, some kind of feedback is required to explain the 
local satellite population. Different power spectra (even 
different values of j3) seem to indicate that different types 
of feedback are needed. For example, in models with 
(Ts^ 0.8 (the precise number depends on typical (3 values 
and the degree of running/tilt), the feedback must be dif- 
ferential. That is, for T^nax — 8 — 30 km s^^, only one out 
of every ~ 5 — 10 halos in this range should form stars. 
On the other hand, in models like the dn/dlnfc = —0.03 
RI model, the BSI case, and our cts — 0.75 model with 
/3 = 0, the discrepancy seems to set in suddenly at Knax ^ 
10 — 20 km , suggesting that nearly all halos smaller 
than this are completely devoid of stars. In this case, the 
feedback mechanism must provide a sharp transition. 

The feedback mechanism proposed by BKW accommo- 
dates the need for only ^ 10% of subsystems to actually 
host observable galaxies by suggesting that only those sys- 
tems that formed before reionization were able to retain 
their gas and eventually form stars. However, if reion- 
ization were to occur very early {e.g., 15), many fewer 
than 10% of these dwarf-sized systems could have collapsed 
before reionization, so that almost all systems smaller than 
Knax ^ 30 km s~^ would be dark. This would be more in 
line with what we sec for the low-power models. This 
is an intriguing result. The best-fit power spectrum of 
the WMAP team (Spergel et al. 2003) leads to simi- 
lar substructure mass fractions as standard CDM, alle- 
viates the central density problem and dwarf satellite dis- 
crepancy, and forces us to consider feedback mechanisms 
that predict a sharp transition between luminous and non- 
luminous galaxies. Additionally, the possible detection 
of early reionization by the WMAP team (Kogut et al. 
2003) provides a feedback mechanism that results in a 
sharp transition. Of course, explaining early reionization 
in models with low small-scale power may be problematic 
(Somerville et al. 2003). Another feedback scenario that 
leads to a sharp transition is photo-evaporation (Barkana 
& Loeb 1999; Shaviv & Dekel 2003). Nonetheless, the un- 
certainty associated with /3 in determining satellite galaxy 
Knax values suggests that efforts to model dwarf galaxy lu- 
minosities as well as dynamical properties will be required 
to resolve this issue (Somerville 2002; Benson et al. 2002). 

Remaining uncertainties in understanding the precise 
nature of the dwarf satellite problem highlight the need 
to focus on attempts to measure CDM substructure by 
other means. Continued efforts to detect siibstructure via 
gravitational lensing {e.g., DKOl; Keeton 2003; Keeton, 
Gaudi, & Fetters 2002; Moustakas & Metcaff 2003) or by 
probes within our own Galaxy {e.g., Johnston et al. 2002; 
Ibata et al. 2002a, 2002b; Mayer et al. 2002; Font et 
al. 2001) offer useful avenues for doing so. Modeling of 
the kind presented here may play an important role in in- 
terpreting results of ongoing observational programs and 
bring us one step closer to confirming or refuting one of 
the fundamental predictions of the CDM paradigm. 
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